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ABSTRACT 

The  problem  of  orbital  parameter  estimation  using  angles  only  observations  is 
examined.  Direction  cosine  measurements,  obtained  from  satellite  passage  of  an 
earth-based  stationary  planar  radar  beam,  are  assimilated  by  an  extended  Kalman 
filter  to  improve  estimates  of  a  classical  orbital  element  set.  Several  progressively 
comprehensive  orbital  motion  models  are  considered  and  compared.  The  relative 
effectiveness  of  these  models  is  illustrated  by  applying  them  to  actual  satellite 
data. 
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I.    INTRODUCTION 

Earth -orbiting  satellites  are  catalogued  in  the  form  of  orbital  element  sets.  For 
maximum  usefulness  of  these  catalogs  to  military  and  civil  space  systems,  element 
sets  need  to  be  updated  as  frequently  and  as  accurately  as  possible.  A  nonlinear 
filter  was  designed  to  perform  these  updates  in  near-real-time  using  observations 
of  satellites  piercing  a  planar  radar  beam. 

The  filter  predicts  time-varying  orbital  parameters  according  to  a  non-linear 
orbital  dynamics  model.  This  model  includes,  along  with  two-body  Keplerian 
motion,  first  order  oblateness  effects,  both  secular  and  periodic,  simple  decreasing 
exponential  atmospheric  drag  effects,  and  white  noise  on  both  the  dynamic  model 
and  the  measured  observations.  This  simple  model  is  found  sufficient  to  track 
several  selected  satellites,  with  proper  tuning  of  the  filter.  It  should  be  noted  that 
even  this  simple  non-linear  model  gives  rise  to  comparatively  complex  filter 
expressions. 

Applying  this  filter  to  the  orbit  improvement  problem  helps  meet  the  need  for 
continual  determination  of  orbit  decay,  collision  avoidance,  satellite  maneuvers, 
etc.  Current  systems  use  a  daily  batch  least  squares  differential  correction  method. 


II.  REVIEW  OF  LITERATURE 

This  research  involves  the  integration  of  a  simple,  yet  adequate,  orbital 
dynamics  model  with  a  necessarily  suboptimal  extended  Kalman  estimator.  This 
combination  is  applied  to  a  geometrically  unique  detection  scheme  to  provide  near 
real-time  improvement  of  orbital  elements  for  selected  satellites.  The  accuracy 
level  under  consideration  falls  under  the  umbrella  of  the  general  perturbations 
problem  of  orbital  mechanics,  which  is  covered  by  numerous  textbooks,  a 
sampling  of  which  are  referred  to  here  [1-6].  In  general,  solution  procedures 
involve  a  gravitational  potential  function  (earth  oblateness,  third  body,  etc.) 
describing  its  contribution  to  the  perturbation  of  a  classical  two-body  Keplerian 
dynamic  model.  The  resultant  expressions  take  the  form  of  variations  in  some  set 
of  orbital  parameters.  The  two  broad  divisions  of  orbital  parameter  representation 
in  the  literature  are:  1)  some  set  of  6  orbital  elements  which  describe  the  orbital 
path  as  a  conic  section  and  its  orientation  in  inertial  space,  and  2)  Cartesian 
position  and  velocity  of  the  orbiting  body.  The  6  classical  orbital  elements  were 
chosen  for  their  physical  significance  in  terms  of  how  the  satellite's  orbit  is  being 
affected  by  various  perturbing  forces.  The  particular  satellite  propagation  theory 
being  used  here  is  loosely  based  on  Brouwer's  theory  [7-11]  for  a  future 
comparison  basis  with  existing  orbit  improvement  procedures.  The  detection 
scheme  is  an  earth-fixed  planar  radar,  the  crossing  of  which  causes  power  to  be 
reflected  back  to  one  or  more  of  multiple  earth-based  receivers  [12]. 


The  most  poorly  modeled  phenomenon  affecting  low  earth  orbits  is  the 
problem  of  atmospheric  drag.  The  overwhelming  consensus  [13-15]  is  that  even 
the  most  complex  atmospheric  density  models  can  only  be  consistently  accurate  to 
within  about  20%.  The  Kalman  filter  incorporates  a  very  simple  decreasing 
exponential  density  model  [16]  corrupted  by  noise  to  track  drag  perturbations. 

The  extended  Kalman  filter  is  used  in  many  situations  for  nonlinear  parameter 
estimation  [17-19]  and  has  been  applied  to  the  field  of  orbit  estimation  in  the  form 
of  orbital  parameter  improvement  by  differential  correction  techniques  [20,21].  It 
is  known  that  such  a  filter  works  well  when  observational  data  are  plentiful 
[22,23].  However,  proper  application  of  such  a  filter  to  a  highly  nonlinear, 
minimal  observation  scenario,  such  as  that  posed  by  the  fixed  radar  plane  detection 
scheme,  can  also  successfully  "track"  satellites. 


III.      TWO-BODY  ORBIT  FORMULATION 

Two-body  orbital  motion  for  a  spherical  earth  can  be  described  by  a  set  of  six 
parameters  (three-dimensional,  second  order  problem).  Traditionally,  these 
parameters  take  one  of  two  forms.  One  method  describes  satellite  motion  in  terms 
of  Cartesian  position  and  velocity.  These  vectors  are  propagated  according  to  the 
two-body  non-linear  equations  of  motion.  The  other  takes  advantage  of  the  fact 
that  the  equations  of  motion  describe  an  orbital  path  which  is  a  conic  section, 
which,  for  earth-orbiting  satellites,  is  always  an  ellipse.  Satellite  motion  can  then 
be  analyzed  as  a  mean  angular  displacement  along  that  elliptical  path.  The  second 
method  is  used  here  because  of  the  linear  nature  of  the  mean  angular  motion  and 
its  clarity  of  physical  significance. 

A.    KEPLER  TWO-BODY  DYNAMICS 

The  orbital  path  is  described  as  an  ellipse  of  fixed  shape  and  size  located  in  a 

plane  of  fixed  orientation  with  respect  to  the  earth's  inertial  reference  frame.   The 

quantities  describing  this  path,  the  classical  orbital  elements  [1],  consist  of  five 

constants  and  one  time-dependent  quantity.   Two  of  the  constants  (a,  e)  describe 
size  and  shape  of  the  orbit,  while  the  other  three  (i,  Q,  (O)  orient  the  orbital  plane  in 

space.    The  time-dependent  quantity  (v)  pinpoints  the  actual  location  of  the 

satellite  on  the  orbital  path.  These  elements  are  (see  fig.  3.1): 
a  (semi-major  axis)  -  defines  size  of  the  elliptical  orbit, 
e  (eccentricity)  -  defines  shape  (non  circularity)  of  orbit. 


Figure  3.1.  Classical  orbital  elements. 


i  (inclination)  -  defines  the  angle  between  equatorial  plane  normal  (K)  and 
orbital  plane  normal  (\V) . 

Q  (longitude  of  the  ascending  node)  -  defines  the  angle  between  the  vernal 
equinox  direction  (I)  and  the  point  where  the  satellite  crosses 
the  equatorial  plane  from  southern  to  northern  hemisphere 
(ascending  node). 

(O  (argument  of  perigee)  -  defines  the  angle  between  the  ascending  node  and 
the  satellite's  closest  point  of  approach  (perigee). 

v  (true  anomaly)  -  defines  the  angle  between  perigee  and  the  satellite's  actual 
position. 

The  elliptical  orbital  path,  is  described  by  the  equation  of  a  conic  section 

a(l-e2) 

r  =  -* -  (3.1) 

1  +  e  cos  v 

where  r  is  the  magnitude  of  the  satellite's  position  vector  at  a  specified  v.    In 
general,  v  is  nonlinear  with  respect  to  tir.ie,  so  the  mean  anomaly  M,  which  varies 

linearly  with  respect  to  time  for  purely  Keplerian  motion,  will  be  used  as  a  filter 
state  instead.  M  is  related  to  v  through  a  quantity  called  eccentric  anomaly  E.  It  is 

defined  as  follows  (see  fig.  3.2): 

M  =  E-esinE  (3.2) 


cosE-e 

cosv  = (3.3) 

1-ecosE 


sinEVl-e2 
sinv  = (3.4) 

1  -  e  cos  E 


Eqs.  3.1-3.4  define  the  relationship  between  M  and  r 
The  filter  state  vector  will  be 


k=[ak     e,     ik     Qt     cok     Mk]T 


(3.5) 


where 


~at+l ' 

ek+; 

= 

wk+1 

Mk+1. 

=  f(xk,Tk) 


(3.6) 


ek 

ik 

Ok 

<*>k 

a~%T  +M. 

i  k  ] 


and  Tk  is  the  time  between  observations.   For  the  classical  2-body  case,  it  can  be 
seen  from  eq.  3.6  that  state  propagation  is  a  linear  problem. 

In  order  to  improve  the  estimate  of  the  state  vector  using  available 
observations,  it  becomes  necessary  to  compare  the  states  with  the  observables  in  a 
common  reference  frame.  Satellite  position,  which  is  implicitly  a  function  of  the 
states  through  the  relationship  in  eqs.  3.1-3.4.  can  be  written 


rp  = 


rcos  v 

r  sin  v 

0 

(3.7) 


where  rp  denotes  a  position  vector  in  the  orbital  (PQW)  reference  frame  (see  fig. 
3.1).  The  remaining  states  are  used  to  transform  satellite  position  from  orbital  to 
inertial  (IJK)  coordinates  using 


Figure  3.2.  Eccentric  vs.  true  anomaly. 


c*  = 


ccocQ  -  scosQci    -scocQ  -  ccosQci      sQsi 
ccosQ  +  scocQci    -scosQ  +  ccocQci    -cQsi 
scosi  ccosi  ci 


(3.8) 


where    c0  =  cos  6    and    sG  =  sin  0  for  brevity.  (See  App.  A  for  development  of 
C».) 


Then  satellite  position  in  inertial  coordinates  will  be 


r1  =C 


X 


rcos  v 

r  sin  v 

0 


(3.9) 


Now,  inertial  satellite  coordinates  need  to  be  related  to  a  measurement  in  the  site- 
based  reference  frame. 


B.    OBSERVATIONS 

The  detection  scheme  modelea  here  is  a  fixed  planar  radar  beam  generated  by 
multiple  continuous  wave  transmitters  along  a  great  circle  path  (see  fig.  3.3).  As  a 
satellite  penetrates  this  beam,  phase  information  is  collected  at  one  or  more  of  the 
receivers.  Antenna  geometry  at  each  site  makes  it  possible  to  compute  direction 
cosines,  which  are  the  pseudo-observations  available  for  element  set  improvement. 
These  observables  are  (see  fig.  3.4) 

p-E 


cos  a  = 


R       PN 

cosp  =  — — 


(3.10) 


where   East-West  cosine  =  cos  a  and  North-South  cosine  =  cos  (3  .    (y  can  be 
calculated  from:   cos2  y  +  cos2  p  +  cos2  a  =  1.) 


We  have 


and 


Figure  3.3.  Radar  plane  geometry 

pcosa  =  p-E 
p  cos  p  =  p  •  N 


z  = 


cos  a 
cosp 


(3.11) 


(3.12) 
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H  (up) 


N  (North) 


->*    E  (East) 


Figure  3.4.  Receiver  site  geometry. 


is  the  vector  of  observables  to  be  used  for  orbital  parameter  improvement.  (See 
App.  B  for  development  of  z  in  terms  of  the  orbital  elements.) 
The  position  vector  in  the  site-based  frame  is 


rH  =CV 


(3.13) 


where 


cH^  = 


c£c§  c£s(J)  s( 

-s<J>c(az)-s£'c<|>s(az)      c(az)c(()  -  s(az)s£s<{>      s(az)c£ 
s(az)s(()-c(az)s^c<j)      -s(az)c<|>-c(az)s£s<t>    c(az)c^ 


(3.14) 


and 


I  =  latitude  of  receiver 
<t>  =  we(t-t0)-lonot£ 
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G)^  =  earth's  rotation  rate 

t0  =  time  of  last  Greenwich  crossing  of  vernal  equinox 

lonjite  =  longitude  of  receiver 

(See  fig.  3.5) 

(See  App.  A  for  development  of  C  " ). 

Recapitulating, 

xk  =f(xk_1,Tk_1) 

zk  =  g(xk,Tk)  (3.15) 

both  of  which  are,  in  general,  nonlinear  functions  of  the  orbital  elements  xk. 
Assuming  for  now  that  perturbations  to  the  two  body  dynamics  and  observation 
noise  can  be  modeled  as  additive  zero  mean,  white  noise  processes,  the  state  and 
measurement  equations  become 

xk  =f(xk.1,Tk_1)  +  wk_1 


(3.16) 


zk  =  g(xk,Tk)  +  vk 

Kalman  filter  equations  are  then 

*kX-i  =f(Xk-X->'T*->) 

pk-,-*p*-k-,*t+« 
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Xw   =  Xk/      +(j,Z,    —  Zw 
7k  7k-l  k  \    k  7*-l 


Pk/k=(l-GkH)P 


« 


k-i 


(3.17) 


Figure  3.5.  Receiver  site  geometry 
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where 


Q  =  E[wkwI] 
R  =  E[vkvJ] 


Pk=Er(xk/{-xk)(xk/k-xk)T 


0  = 


3f(xk,Tk) 


3x, 


'X-i'1"'- 


H  = 


^g(xk,Tk)| 


ax 


k         ,*X_,-T«- 


£x-issfi(x«;-i'Tw) 


zk  = 


pcosa 
pcosp 


Jk 


(3.18) 


(See  App.  B  for  derivations  of  <P  and  H.) 

In   words,   the  actual   observations,    zk,   are   compared   with   computed 
observations,  iy     (based  on  the  classical  two-body  propagation  model),  resulting 

in  the  Kalman  filter  innovations.  These  are  multiplied  by  the  filter  gain,  which  is 
based  on  the  dynamic  model  and  measurement  uncertainties  and  used  to  produce 
an  improved  estimate  of  the  orbital  elements. 
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C.    RADAR  PLANE  INTERSECTION 

Meaningful  application  of  the  observations  to  orbital  element  set  improvement 
requires  the  ability  to  predict  subsequent  penetrations  of  the  radar  plane  by  the 
satellite  in  question.  A  simple  way  to  do  this  is  to  define  the  satellite's  position 
vector  in  a  coordinate  frame  referenced  to  the  radar  plane  (See  fig.  3.6).  Then  the 
out  of  plane  component  (Z)  of  position  can  be  checked  for  a  value  of  zero, 
yielding  an  in-plane  condition. 

Satellite  position,  readily  available  in  orbital  frame  coordinates  (rp)  is 
transformed  to  radar  plane  coordinates  by  performing  two  coordinate 
transformations:  one  from  orbital  to  inertial  frame,  the  other  from  inertial  to  radar 
plane  frame,  or 


rx  =Cx^C^rp 


(3.19) 


where  rx  denotes  a  position  vector  in  the  radar  plane  (XYZ)  reference  frame  and 


c^  = 


c8  ce 

-s6 


cese    se^ 

rp  rp 


c9 


0 


-s0   c6     -s0    s9    c9 

rp  rp  rp 


(3.20) 


where 

e  =  coe(t-t0)-lonx, 

8^  is  the  inclination  of  the  radar  plane  to  the  equatorial  plane 
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Ion,  is  the  longitude  where  the  maximum  latitude  of  the  radar  plane  occurs, 
which  describes  the  location  of  the  X-axis  of  the  radar  plane 
reference  frame.  (See  fig.  3.6). 

(See  App.  A  for  development  of  C  ^ .) 


Figure  3.6.  Radar  plane  orientation. 

Since  the  only  component  of  satellite  position  of  interest  in  determining  radar 
plane  intersection  is  the  Z-component,  the  expression  can  be  further  simplified  to 
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rz=[0     0     l]Cx^C^rp 

=  [-seipP1c6-seipP4s0  +  cerpP7]rcosv 

+[-86^06 -s0ipP5s6  +  ceipPl]r  sin  v 


(3.21) 


where 


Q%   = 


P,        P, 


P         P 

r  7         r  8 


9  _ 


(3.22) 


The  in-plane  condition  of  relevance  to  the  filter  is  that  which  falls  inside  the  field 
of  view  of  one  of  multiple  receiver  sites,  also  placed  along  the  great  circle  path 
circumscribed  by  the  transmitters.  The  time  epoch  of  this  condition  is  the  time 
close  to  which  we  expect  to  observe  the  satellite  from  one  or  more  of  the  radar 
plane  receivers. 
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IV.      DOMINANT  PERTURBATIONS 

Satellite  orbital  parameters  deviate  from  Keplerian  motion  deterministically 
and  randomly  due  to  many  factors.  Ignoring  for  now  such  intentional  human- 
caused  phenomenon  as  plane  change  and  stationkeeping  maneuvers,  the  largest 
natural  contributor  to  parameter  variations  is  the  nonuniformity  of  the  gravitational 
field  due  to  earth's  oblateness.  For  low  earth  orbits,  the  next  most  significant 
factor  is  atmospheric  drag. 

A.    GRAVITY  PERTURBATIONS 

The  largest  contributor  to  perturbations  to  the  Keplerian  orbital  elements 
arises  from  the  oblateness  of  the  earth.  These  perturbations  take  the  form  of  a 
combination  of  secular  and  periodic  variations  of  the  classical  orbital  elements. 
Many  solutions  have  been  developed  for  this  problem,  each  with  its  own  strengths 
and  weaknesses.  All  the  solutions,  however,  are  based  on  the  gravitational 
potential 

U  =  ^{l-i^J„P„(sin5))  (4.1) 

where 

fie  is  earth's  gravitational  parameter 

r    is  the  distance  of  the  satellite  from  earth's  center 

Jn  are  n^  order  zonal  harmonics 
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Pn(sin  5)  are  Legendre  polynomials 
5  is  satellite  declination 

and  are  more  or  less  complex  depending  on  the  order  to  which  the  terms  of  the 

expression  are  carried.   This  expression  for  U  already  carries  with  it  a  degree  of 

simplification  in  that  it  assumes  axial  symmetry  for  the  earth.  In  truth,  the  equator 

is  slightly  elliptical,  but  the  main  effect  is  on  the  stability  of  geosynchronous 

satellites,  since  their  orbits  are  in  the  equatorial  plane  [3].  Further  simplification  is 
accomplished  by  dropping  terms  for  n  >  3.   This  is  justified  for  medium  accuracy 

applications  because  of  relative  magnitudes  of  the  Jr  terms 

J,  =  1.0826 -10"3 


J3  =-210 


J5  =-210 


-t. 


-7 


(4.2) 


where  Jn  successively  decrease.  So  the  simplified  potential  function  is 

u  =  _fH±  +  ^j2(3sin26_1)'|  (43) 

Defining  R  =  — f-  J 2 (3 sin2  6-l)  as  the  perturbing  function,  the  variation  of  that 
function  with  respect  to  each  of  the  orbital  elements  is  [  4  ] 
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dR  _ 

aa " 

-1r 

a 

dR 

3J2 

r         /■■ 

•)(!• 

ae  ~ 

2(l-e2)r 

-lcosco(l  +  ecos,\ 

-sin2  i  sin 

2usin  v(2  +  ecos  v)] 

dR  _ 

-3J2 

.         .       "> 

di  " 

— - — sin  i 

•cosi  sin^  u 

aR 

-3J2       2 

aco 

—  sin 

2r3 

i  sin  2u 

aR 

an" 

0 

aR 

1  dR 

8M 

n  dt 

(4.4) 


To  develop  a  first  order  theory,  the  right-hand  members  of  eq.  4.4  are 

expanded  in  a  Taylor  series  about  the  initial  orbital  element  values,  retaining  only 

the  first  term  of  the  expansion,  yielding  expressions  for  the  derivative  of  each 
orbital  element  with  respect  to  true  anomaly  v.    Integrating  these  expressions 

yields  an  analytical  first  order  solution  in  the  form 

x  =  xo+5sx(v)  +  6px(v)  (4.5) 

where 

x  =  instantaneous  orbital  elements 

x0  =  initial  orbital  elements 

8sx  =  secular  variations  as  a  function  of  v 

5pX  =  periodic  variations  as  a  function  of  v 
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All  six  orbital  elements  experience  periodic  variations.    However,  only  Q,  co 
and  M  experience  secular  perturbations.  These  are 

5Q  =  -^t[cosio]2(v-vo) 
4p 


s  3J? 

5,co  =  — r 


2p" 


2 —  sin^  i 

9 


(v-vo) 


(4.6) 


8  M  =  n 


J,  (\ 


V 


i3A 


1  +  -T- 

aolro; 


(l-3sin250) 


(t-tj 


Note  that  6SM  is  expressed  in  terms  of  independent  variable  t,  which  is  the  time  at 
which  true  anomaly  is  equal  to  the  value  of  v  in  5sQ  and  §sco.  The  periodic 
variations  in  co  and  M  do  cam'  l/e0  terms  which  approach  singularity  for  near 
circular  orbits.   This  can  be  explained  by  the  fact  that  co  and  M  are  undefined  for 

circular  orbits  and  ill-defined  for  near  circular  orbits.  Fortunately,  the  combination 
of  co  +  M  obviates  the  singular  condition  by  effectively  forming  one  new  orbital 

element  out  of  co  and  M. 

Long  term  propagation  of  the  orbital  elements,  assuming  a  drag-free 
environment,  is  accomplished  by  applying  only  the  secular  variations,  since 
periodic  variations  will  have  no  net  effect  on  the  elements  over  time.  However, 
the  periodic  variations  will  be  seen  later  to  come  into  play  when  observations  are 
assimilated  into  the  orbital  element  improvement  process.  Inclusion  of  these 
secular  variations  in  the  dynamic  model  yields 


21 


xk+,  =  f(*k'TJ 


(4.7) 


where 


f(xk,Tk)  = 


3J, 


Q*  -T-T^OOK^i-vJ 


2p: 


co.  + 


3J„ 


^     2Pk2 


2  —  sin'  i, 

i  k 


(vk+1-vk) 


Mk  +  nk 


J2 

1+? 


^o        > 


vr.  y 


-3/ 


(l-3sin26k) 


T„ 


(4.8) 


where 

Pk  =  ak  (l  — ek )  is  the  semi-latus  rectum  of  the  orbit 

J2  is  the  second  order  zonal  harmonic  (equatorial  bulge). 

The  gravitational  force  model  which  gives  rise  to  these  secular  variations  is 
relatively  simple.  It  assumes  that  the  earth  is  symmetric  about  its  spin  axis  and 
about  the  equatorial  plane.  Noteworthily,  the  first  three  mean  orbital  elements 
experience  no  secular  variation  due  to  earth's  oblateness. 
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It  may  be  possible  to  maintain  identification  of  some  satellites  using  secular 
variations  alone.  Continuous  identification  of  a  satellite,  however,  may  depend 
upon  the  application  of  the  periodic  variations  to  the  mean  elements  at  observation 
times.  These  periodic  variations  are  assimilated  as  an  additive  (superposition 
assumed)  adjunct  to  the  existing  dynamic  model  and  therefore  are  not  technically  a 
part  of  the  propagation  dynamics.  The  filter  innovations  resulting  from  a 
comparison  of  observables  with  predicted  states  are  then  applied,  via  the  Kalman 
gains,  to  the  mean  elements,  implying  a  prior  removal  of  the  periodic  variations. 
These  periodic  variations, to  first  order,  and  ignoring  terms  of  the  order  of  e0  and 
smaller,  are 

3J.,a. 


=  -2»o 

P      2p02(l-e02) 


sin:  i0-cos[2(co0+v)] 


3J2 
oe    = 

r     2p;L 

7 


1  — sin2  io    cos  v  +  —  sin2  io  -cos(2cdo  +  v) 


+  —  sin2  io  -cos(2gl>o  +3v) 


8i    —Hi 


3  J 


sin(2io)cos[2(coo+v) 
2(co0  +  v) 


cosi     sin 


5Qp  =— ^cosio  •sin[2(coo  +  v)] 

Mr  o 


(4.9) 


3J    ( 

P         4p2l 


1--sin2ioJsin[2(coo+v)] 
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B.    DRAG  PERTURBATIONS 

Satellites  traveling  in  the  altitude  regime  below  600  km  experience  non- 
negligible  drag  force  due  to  atmospheric  resistance  [3].  This  force  is  well 
represented  as 

FD=^CD-pV2  (4.10) 

2        m 

where 

m=mass  of  satellite 

CD~2  is  a  dimensionless  drag  coefficient 

A=  effective  cross  sectional  area  of  satellite 

V=  velocity  of  satellite  with  respect  to  atmosphere 

p=  local  atmospheric  density. 

1.      Atmospheric  Drag  Effects 

To  first  order,  drag  force  doesn't  affect  i  and  Q  [6].   The  main  effects  of 

this  drag  force  are  to  bring  about  secular  decreases  in  a  and  e,  i.e.,  to  decrease  the 
size  and  eccentricity  of  the  orbit,   oo  and  M  also  experience  changes  due  to  drag, 

but  these  are  negligible  compared  to  the  oblateness  effects. 
It  turns  out  that 

da_     CDApaVl  +  ecosE^(1_T)J  (4n) 


dM  m  Vl-ecosE. 


where 
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d(l  -ecosE) 

t  = 

1  +  ecosE 


n    <  '     ^  (4.12) 


d  =  ^L(!_e^ 


COSl 


V- 


and 


de  CDA      ,  (1  +  ecosE)^  2 

= — pa(l-e   ) 37(1-1)    cosE 

dM  m  (l-ecosE)/2 

(4.13) 

CDA        /1-ecosEV*         .   .   2_ 

— paed  (l-x)sin    E 

2m  U  +  ecosEy 


All  quantities  in  these  two  expressions  are  well  known  (according  to  the  simplified 
model)  except  for  atmospheric  density  p,  the  value  of  which  is  related  to  space  and 

time  in  a  complex  manner.    Atmospheric  models  delivering  as  good  as  15% 

standard  deviation  are  in  the  form  of  tables  or  very  complicated  empirical 

formulae  or  both  [13].    For  simplicity,  it  is  advantageous  to  use  an  analytical 
function  to  model  p.    It  is  assumed  that  variations  in  p  in  the  upper  atmosphere 

(150  km  to  750  km)  are  a  superposition  of  four  individual  variations  [4]. 

The  first  variation  is  diurnal,  or,  daily.  This  fluctuation  has  to  do  with 
change  in  solar  flux  intensity  associated  with  moving  from  darkness  to  daylight, 
and  vice  versa.  Density,  depending  on  latitude  and  altitude,  may  change  by  up  to  a 
factor  of  10. 
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The  second  type  is  related  to  the  11  year  solar  cycle  which  is  motivated 
by  intensity  of  sunspot  activity.  Maximum  density,  occurring  at  the  peak  of 
sunspot  activity,  is  roughly  ten  times  the  mean  density. 

A  third  variation  is  semiannual,  with  density  falling  to  a  minimum  in  July 
and  a  maximum  in  October,  with  a  less  pronounced  minimum  and  maximum  in 
January  and  April,  respectively.  The  range  of  this  variation  is,  generally,  less  than 
1/20  the  range  of  the  1 1  year  solar  cycle. 

The  fourth  type  of  variation  is  nonperiodic  and  unpredictable  and  is  tied 
to  geomagnetic  storms  rising  from  solar  flare  eruptions.  Atmospheric  density  may 
increase  by  a  factor  of  10  from  the  geomagnetic  quiescent  value. 

Over  the  short  term  (days),  barring  magnetic  storms,  the  diurnal  variation 
dominates.  However,  over  several  years,  the  solar  cycle  has  a  more  profound 
effect  on  atmospheric  density. 

The  modeling  approach  followed  here  was  to  start  with  a  very  simple, 
even  crude,  atmospheric  density  model,  including  the  above  considerations  as 
necessary  to  achieve  desired  accuracy. 

2.      Simplified  Drag  Effects 

First,  simpler  expressions  for  changes  in  a  and  e  due  to  drag  can  be 
developed  by  viewing  these  changes  as  a  result  of  orbital  specific  energy  being 
extracted  by  drag  force  [14].  Specific  energy  of  an  orbit,  which  is  conserved  for 
an  unperturbed  Keplerian  orbit,  can  be  expressed  as 

£  =  i-V2-^  =  -  —  (4.14) 

2  r         2a 
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The  change  in  a  for  an  incremental  change  in  s  is  then 


2a 
da  = de 


C4.15) 


If  it  is  assumed  that  most  of  the  drag  occurs  near  perigee  (atmospheric  density 
exponentially  decreasing  with  altitude),  then 


r  =  r   =a(l-e) 


(4.16) 


from  eq.  3.1,  and 


v     ay 
=  Ji(l  +  e) 


(4.17) 


from  eqs.  4.13  and  4.15.   For  a  satellite  moving  through  dv  near  perigee  (see  fig. 
4.1),  the  specific  energy  extracted  will  be 


Figure  4.1.  Satellite  moving  through  perigee. 
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dE  =  -5Lrdv 

m 

= 5_  pv2rdv 

21    m    j 


(4.18) 


Substituting  eq.  4.17  into  eq.  4.18  yields 


de-I 

2 


rc.A^ 


V  m  ; 


pfi(l  +  e)dv 


(4.19) 


Substituting  eq.  4.19  into  eq.  4.15  gives  an  expression  for  decrease  in  semi-major 
a^is  due  to  a  decrease  in  orbital  energy 


da  =  -|-D 

V    m    j 


a:(l  +  e)pdv 


(4.20) 


Since  the  new  orbit  will  pass  through  the  same  perigee  point,  all  of  da  shows  up  as 
a  decrease  in  apogee  radius,  and 


dr   =  da(l  -  e)  -  ade  =  0 


(4.21) 


by  differentiating  eq.  4.16.  Solving  for  de  gives 


a 


— ^—  a  1-e2  pdv 
V    m    J 


(4.22) 


3.      Atmospheric  Drag  Characterization 

A  crude  model  for  variation  of  p  starts  with  the  assumption  that 

temperature  and  chemical  composition  of  a  gas  remain  constant  [14].    Then, 
density  is  directly  proportional  to  pressure  according  to 
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p  =  kP  (4.23) 

(Realistically,  k  is  a  function  of  temperature  and  chemical  composition.)  At  an 
altitude  h,  the  decrease  in  pressure  accompanying  an  increase  in  altitude  is 

dP  =  -pgdh  (4.24) 

where  g  is  acceleration  of  gravity.  Since  k  is  assumed  constant. 

dP  =  -dp  (4.25) 


Combining  eqs  4.24  and  4.25  yields 


dp 

-^  =  -kedh 


dp 


J— =  -kgfdh  (4.26) 

P.    P  o 


P  =  Poe'k*h 


So,  in  general,  atmospheric  density  decreases  exponentially  with  altitude. 
Eq.  4.26  can  be  rewritten 

p  =  poe"X.  (4.27) 

where  p0  (reference  air  density)  and  h0  (density  scale  height)  change  for  different 
altitude  regimes.   For  example,  using  standard  fixed  values  of  p  from  the  "U.S. 
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Standard  Atmosphere,  1976,"  a  reasonable  fit  can  be  made  by  applying  the  values 
in  Table  4.1  for  the  ranges  of  altitude  shown. 

Table  4.1.  DRAG  PARAMETERS. 


Altitude 

Po(kg/m?) 

z0  (km) 

80km<z<120km 

18.5 

5.8 

120km<z<140km 

1.21x  10"3 

11 

140km<z<200km 

3.00  x  lO"6 

21 
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V.   THE  RECURSIVE  ORBITAL  PARAMETER  ESTIMATOR 

An  extended  Kalman  filter,  designed  for  orbital  parameter  estimation,  is  based 
on  the  earth-fixed  radar  plane  detection  scenario  described  in  Chap.  Ill  (Also  see 
fig.  5.1).  The  filter  propagates  the  orbital  elements  in  time  according  to  a  first 
order  nonlinear  dynamic  model  which  takes  into  account  the  dominant 
perturbations  arising  from  earth  oblateness  and  atmospheric  drag.  Upon  receiving 
observations,  in  the  form  of  pairs  of  direction  cosines,  the  filter  calculates  an 
improvement  to  the  current  parameter  estimate.  The  filter  operates  recursively  as 
observations  become  available.  The  filter  is  written  in  Matlab  code  and  can  be 
found  in  App.  C. 

A.    FILTER  INITIALIZATION 

The  estimator  is  designed  to  receive  two  data  files.  One  contains  an  initial  set 
of  orbital  elements  along  with  subsequent  pairs  of  available  direction  cosines  (see 
App.  C).  The  other  holds  information  about  the  receiver  sites.  Values  of  pertinent 
constants  are  also  set,  as  are  the  initial  times  necessary  for  orbit  propagation  and 
estimation  of  observables. 
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INPUT  DATA 


INITIALIZE  KALMAN  FILTER 


PROPAGATE  STATES 


ITERATE  TO  IN-PLANE  CONDITION 


ESTIMATE  MEASUREMENTS 


CALCULATE  <t>,P,Q,H,G 


APPLY  STATE  IMPROVEMENTS 


YES 


Figure  5.1.  Estimator  flow  chart. 
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1.      Input  Data  Files 

Two  input  files  are  required  for  filter  operation.  One  contains  an  initial 
set  of  orbital  parameters,  along  with  subsequent  sets  of  direction  cosines.  This 
information  is  all  tagged  with  corresponding  times  (see  App.  C).  The  last  five 
orbital  elements  occur  in  classical  orbital  element  form,  but,  in  lieu  of  semi-major 
axis,  orbital  period  is  given.  The  conversion  is  simply 


a  = 


1^ 

271 


[in  J 


(5,1) 


where  Tsalis  the  period  of  the  satellite.   Note  here  the  apparent  disappearance  of 
earth's  gravitational  parameter  }ie.  This  convenience  is  brought  about  by  the  use 

of  canonical  units  [1],  where 


H*  =1 


TUi 


1  DUe=  6378.135  km  and  is  earth's  mean  equatorial  radius. 

1   TUe=  13.44686  min  and  is  the  time  it  takes  for  a  satellite  traveling  in  a 
circular  orbit  at  a  radius  of  1  DU  to  travel  1  DU  of  arc  length. 

Canonical  units  are  used  in  all  filter  calculations.    The  other  input  file  contains 

information  on  the  receiver  site  locations, altitudes,  and  azimuths.   The  azimuthal 

measurement  indicates,  for  each  receiver,  the  angular  offset  of  its  coordinate  frame 

from  true  north. 
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2.      Constants 

The  earth's  apparent  rotation  rate  is  initialized  as 

coe  =.05883378    ra%u  (5.2) 

The  radar  plane's  inclination  to  the  equatorial  plane  and  the  longitude  of  its  center 
line  are,  respectively, 

0^  =33.58310° 
Ion.  =-101.31348° 


(5.3) 


The  radar  plane  is  offset  from  earth's  center  by  about  20  km.  or 

roflsel=.0031    DUe  (5.4) 


3.      Time 

All  times  are  given  in  Universal  Coor  jinated  Time  (UCT)  and  converted 
to  TU©s  for  filter  use.  Tgn^h. which  is  time  elapsed  since  the  Greenwich  meridian 
last  passed  the  geocentric  inertial  I  axis,  is  computed  from  a  reference  angle 
(measured  from  I  to  Greenwich)  at  0  hrs,  1  Jan  92,  of 


0.  =  98.92025093  lc  (5.5) 


using  the  coe  given  previously 
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4.      Filter  Matrices 

The  error  covariance  matrix  P  is  initialized  to  reflect  the  level  of 
confidence  placed  in  the  accuracy  of  the  initial  orbital  element  set.  The 
measurement  error  covariance  matrix  R  is  set  to 


R  = 


cosa 

0       a 


0 

2 


(5.6) 


where  oco6a  =ocotP  =.0002°  reflects  the  uncertainty  of  the  available  direction 
cosines. 


B.    STATE  PROPAGATION 

The  filter  states  are  propagated  according  to  the  two-body  Keplerian  orbital 
dynamics  perturbed  by  earth  oblateness  and  atmospheric  drag.  The  nonlinear  state 
equations  are 

CDAla2(l  +  e)p(vk+:-vk) 


f(*k,Tk)  = 


V    m    j 


ek  " 


n,  - 


m   ; 
3J 


2pk 


a(l-e2)p(vk+{-vk) 
[cos(ik)](vk+1  -vk) 


CO,  + 


3J 


k     2pk2L 


5.2. 

2  —  sin    i. 

2 


(vk+1-vk) 


Mk  +nk 


^ 


h(*A 


-,3/ 


vrw 


(l-3sin26k) 


(5.7) 
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and  are  propagated  at  arbitrary  increments  of     "yg  .    This  is  continued  until  a 

crossing  of  the  radar  plane  is  detected  by  a  sign  change  of  the  radar-plane-normal 
component  of  satellite  position. 

C.    RADAR  PLANE  INTERSECTION 

Detection  of  a  radar  plane  crossing  calls  a  function  which  iterates  to  an  in- 
plane  condition  for  estimated  satellite  position.  This  is  accomplished  using  a 
Newton-Rapson  iteration  on  the  expression  for  the  plane-normal  component  of 
satellite  position.  The  iteration  calls  for  a  calculation  of  an  approximate  AM 
which  will  drive  rz  to  near  zero.  This  is  given  by 

AM=      ~r'    N  (5.8) 

/dM 


where,  from  eq.  3.21, 


r   =  X.rcos  v  +  X,rsin  v  (5.9) 


and 

dr       dX,  3(rcosv) 

— -  = Lrcosv  +  X,  — 

BM      BM  dM 

9X,      .  3(rsinv) 

+ rsmv  +  X, 

dM  2       dM 

Taking  the  individual  derivatives  yields 


(5.10) 
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where 


and 


where 


ax.         .  „   a  cose    „   .  n  aisme 

L  =  -P.  sin  6    — -  P4  sin  e„  — 

a\l  l         *      dM  4         '     8M 


ax2  a(cose)       .      a(sme) 

=  -P,  sin  8m Pc  sinG 

dM  *      dM  5  *     3M 


a(cose)         _3/, .  _, 

=  -coea  ^(sin6) 
dM 

a(sin6)  _y 

— — —  =  CDe2    ^(COS0) 

a\i 


a(rcosv)  asinE 


a.M  1  -  e  cos  E 


a(rsin  v)         a\l  -e:  cosE 
dM  1-ecosE 


(5.11) 


(5.12) 


(5.13) 


rcos  v  =  a(cosE-e) 

rsin  v  =  aVl-e:  sinE  (5.14) 

aE  1 


dM      1-ecosE 

This  iteration  is  found  to  converge  to  within  zr  <e  in  three  to  four  iterations, 
where  e  =  10~8. 
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The  function  then  determines  whether  the  satellite  is  within  the  field  of  view 
of  one  or  more  receiver  sites  by  calculating  the  in-plane  rx  and  ry  components  of 
satellite  position  and  checking  the  inequality  (see  fig.  5.2) 

r. 


F: 


> cos  20' 


(5.15) 


V      x  > 


It  should  be  noted  here  that  the  field  of  view  needed  to  be  opened  up  to  pick  up 
some  satellite  observations  by  receivers  at  the  edges  of  the  radar  plane. 


TX  (Ion  =101°) 


Y 


Figure  5.2.  Sensor  field  of  view. 

D.    ASSIMILATION  OF  OBSERVATIONS 

Once  a  satellite  is  in  the  field  of  view  of  radar  plane  coverage,  the  input  data 
file  is  checked  for  observations  at  or  near  its  time  of  arrival.  If  a  pair  of  direction 
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cosines  is  available,  then  a  measurement  estimate  is  calculated  according  to  the 
development  in  App.  B,  where 


z  = 


Pn 


(5.16) 


O  is  also  calculated  according  to  App.  B.  The  plant  noise  covariance  matrix  Q  is 
calculated  using  the  squares  of  expected  levels  of  uncertainty  in  each  of  the  states. 
For  example,  if  the  two-body  problem  is  the  dynamic  model,  then  the  Q  matrix 
would  be 


Q  = 


?.2 

0 

0 

0 

0 

0 

0 

9 
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0 

0 

0 

0 

0 

G2 

1 

0 

0 

0 

0 

0 

0 

< 

0 

0 

0 

0 

0 

0 

< 

0 

0 

0 

0 

0 

0 

G2 

(5.17) 


where  Gnare  proportional  to  the  magnitude  of  the  largest  unmodeled  perturbation, 
in  the  two-body  case,  earth  oblateness.  So,  for  the  two-body  problem, 
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r°' 

°< 

Oi 

°Q 

<>■ 

L^mJ 

3J,a 


(l  -  e 2 )  sin :  i  •  cos[2(co  +  v)] 


3J,  r.  3  . 


^P' 


1 


1 — sirT  i   cos  v  +  — sin*  i-cos[2(o  + vl 

V       2  J  4  L  J 

3J 

— 2-  sin  2i  •  cosf2(co  +  v)l 
8p2  L  J 

3J, 

-  n  cos  i  •  T, 


2p 
3J,     ( 

2p2 


2  —  sin2  i  T, 


V 


— Tn  l--sin    i    1-e        T 
2p2     I      2  /  ' 


(5.18) 


Now,  the  error  covariance  matrix  can  be  calculated  as 


Pk/     =ip        <fcT+Q 


A-\ 


k-l 


(5.19) 


H  is  then  calculated  as  given  in  App.  B,  making  possible  the  calculation  of  filter 
gain 


Gk=P.    FT(HP.    HT+R 


A-l 


(5.20) 


which  is  then  applied  to  the  residuals   z-z  to  calculate  improved  orbital 
parameters  \y. 

It  should  be  noted  that,  throughout  the  filter,  the  quantities  v  and  E  are  needed 
for  various  calculations.  E  is  calculated  fom  M  using  Newton-Rapson  iteration  on 


M  =  E  -  e  sin  E 


(5.21) 


and  v  is  calculated  by  solving  for 
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v  =  sin 


-i 


VI- e2  sinE 
1  -ecosE 


(5.22) 


then  using 


V  =  cos 


"  cosE-e 
_l-ecosE. 


(5.23) 


to  resolve  the  ambiguity 
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VI.      RESULTS 

This  extended  Kalman  filter  is  designed  to  estimate  the  orbital  parameters  of  a 
satellite,  given  observables  in  the  form  of  a  pair  of  direction  cosines  obtained  upon 
the  satellite's  passage  of  an  earth-fixed,  stationary  planar  radar  beam.  One  day's 
worth  of  data  for  each  of  three  satellites  was  obtained  for  filter  testing  and 
adjustment.  Results  of  the  application  of  three  successively  more  comprehensive 
filter  models  are  compared  with  each  other. 

The  first  filter  models  only  two-body  motion,  the  second  adds  secular 
perturbations  due  to  earth  oblateness.  and  the  third  also  includes  periodic 
perturbations  due  to  oblateness  as  well  as  atmospheric  drag  effects.  Each  filter,  of 
course,  models  unknown  perturbations  as  zero  mean,  white,  Gaussian  noise,  in  the 
tradition  of  Kalman  filtering.  Unfortunately,  none  of  the  available  satellite  data 
include  an  orbit  in  the  drag  regime.  The  filters  will  be  compared  in  the  areas  of 
radar  plane  time  of  arrival,  a  residual  small  sample  statistic,  and  orbital  parameter 
estimation. 

A.    RADAR  PLANE  TIME  OF  ARRIVAL 

With  orbiting  objects  arriving  at  the  radar  plane  on  the  order  of  every  few 
seconds,  it  is  important  that  prediction  of  arrival  times  be  fairly  accurate.  Current 
methods  yield  accuracies  of  within  a  second.  The  best  that  can  be  hoped  for  under 
the  current  detection  configuration  is  about  0.25  second,  which  is  the  sensor  time 
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uncertainty.   After  some  crude  tuning  of  the  three  filters,  a  comparison  of  arrival 
times  for  the  three  satellites  is  shown  in  Table  6.1,  where  At  =  t  .    -t 


ok  es: 


It  can  be  seen  that,  in  all  cases  and  for  all  models,  the  first  detection  time  after 
the  estimator  is  started  up  (which  occurs  on  the  order  of  12  hours  after  filter 
initiation)  is  always  mispredicted  by  the  filter  by  a  time  on  the  order  of  seconds. 
However,  it  is  seen  that  the  secular  and  periodic  models  cut  down  that  time  error 
by  a  factor  of  from  three  to  four.  Subsequent  detections  yield  differences  in  time 
between  calculated  and  observed  plane  intersections  of  less  than  1  second  for  a 
one-orbit-later  detection  and  between  two  and  three  seconds  if  detection  occurs 
next  on  the  "backside"  of  the  orbit  (about  a  half  day  later).  Although  this  is  true 

Table  6.1.  OBSERVED  VS.  ESTIMATED  TIMES  OF  ARRIVAL. 


2-bodv 

Secular 

Periodic 

At  (sec) 

Satellite  116 

13.94 
2.21 

4.29 
0.23 

4.30 
0.58 

-3.64 

-2.61 

-0.63 

-0.32 

0.32 

0.87 

Satellite  117 

15.18 
0.78 

4.76 
-0.10 

4.78 
0.17 

-0.13 

-0.18 

0.27 

-0.69 

-2.40 

-0.93 

Satellite  118 

12.39 
0.63 

3.49 
-0.01 

3.49 
0.21 

-2.24 

-2.03 

-0.81 

0.32 

0.25 

0.59 
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for  the  two-body  and  secular  cases,  the  incorporation  of  periodic  perturbations 
yields  subsequent  time  differences  of  less  than  one  second  in  all  cases. 

B.    RESIDUAL  SMALL  SAMPLE  STATISTIC 

The  filter  residual  is  the  difference  between  actual  and  estimated  observations, 
and,  as  such,  is  the  basis  for  a  good  performance  measure.  The  small  sample 
variance  of  the  residual  is  calculated  for  each  case  as  follows: 

G2=£--X2  (6.1) 

.=:    n 

where 

X  =  X—  (6.2) 

i-i   n 

is  the  small  sample  mean,  x,  being  the  i^  residual.  Plots  of  these  variances  for  each 
case  are  shown  in  figs.  6.1  -  6.3. 

It  is  evident  that  the  residuals  for  the  N-S  direction  cosines  are  always  very 
small.  This  is  because  this  direction  cosine  is  always  close  to  90°  and  the 
geometry  of  the  problem  always  forces  the  estimate  of  this  direction  cosine  to  be 
very  close  to  90°.  However,  the  E-W  (in-plane)  cosines  are  more  interesting  to 
observe.  In  general,  the  trend  is  toward  smaller  residual  variances  from 
observation  to  observation  and  from  simple  to  more  complex  dynamic  model.  The 
only  case  where  this  does  not  hold  true  is  for  Satellite  #117,  where  it  can  be  seen 
that  the  two-body  model  is  sufficiently  lost  after  approximately  eight  hours  to 
yield  higher  residuals  than  at  the  previous  set  of  observations.  Note  that  a  "stack" 
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of  asterisks  denotes  several  observations  gathered  at  the  same  radar  plane  crossing, 
which  are  then  applied  one  at  a  time  for  parameter  improvement. 

It  can  be  concluded  from  this  limited  residual  data  that,  for  the  relatively  short 
times  included  in  the  data,  the  filter  is  performing  acceptably. 

C.    ORBITAL  PARAMETER  ESTIMATION 

The  ultimate  objective  of  this  filter  is  to  be  able  to  predict  what  a  satellite's 
position  will  be  at  some  future  time.  Hopefully,  the  filter  is  predicting  plane 
intersection  times  by  accurately  estimating  orbital  parameters.  Figs.  6.4  -  6.6 
show  time  propagating  values  of  the  four  parameters  that  are  expected  to  be 
varying  most  with  respect  to  their  two-body  propagation  values.  All  three 
satellites  show  similar  results.  In  all  cases,  semi-major  axis  estimates  differ 
greatly  from  model  to  model.  Eccentricity  estimates  for  the  two-body  and  secular 
models  are  very  similar,  with  the  periodic  model  showing  the  greatest  deviation. 
All  estimates  are  very  similar  for  the  longitude  of  the  ascending  node,  though, 
between  observations,  the  two-body  model  does  not  have  a  mechanism  to 
accurately  propagate  this  element.  Argument  of  perigee  estimates  are  very  close 
between  the  secular  and  periodic  models,  while  the  two-body  makes  very  little 
correction  to  this  parameter. 
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Figure  6.1.  Variance  of  residuals  (sat.  #  116). 
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Figure  6.2.  Variance  of  residuals  (sat.  #  117). 
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Figure  6.3.  Variance  of  residuals  (sat.  #  118). 
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Figure  6.4.  Orbital  element  estimates  (sat.  #  116). 
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Figure  6.5.  Orbital  element  estimates  (sat.  #  117). 
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Figure  6.6.  Orbital  element  estimates  (sat.  #  118). 
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VII.  CONCLUSIONS  AND  RECOMMENDATIONS 

Orbital  parameter  estimation  was  attempted  using  an  extended  Kalman  filter 
for  each  of  three  orbital  dynamic  models.  Observations  were  angles  only  and 
made  roughly  twice  a  day  per  satellite  passing  through  the  radar  plane.  It  was 
found  that  the  two-body  orbital  model  including  secular  and  periodic  variations 
due  to  earth  oblateness,  was  able  to  predict  subsequent  radar  plane  crossings  to 
within  a  second  of  actual  crossings. 

The  filter,  which  includes  a  crude  model  for  atmospheric  drag  effects,  was  not 
tested  on  any  satellites  in  the  drag  regime.  Future  research  should  include  data 
from  satellites  with  altitudes  less  than  600  km.  Also,  longer  term  data  should  be 
obtained  to  check  filter  performance  over  the  long  haul  and  also  to  see  if  the  filter 
could  track  longer  term  disturbances  or  perturbations  to  the  orbit  under 
consideration. 

Tuning  of  the  filter  could  also  be  accomplished  by  obtaining  a  larger  data  base 
of  satellite  orbital  elements  and  corresponding  observations.  This  could  be  done 
manually,  but  it  would  be  more  exciting  to  see  if  an  artificial  neural  network  could 
be  designed  and  trained  to  accomplish  the  same  task. 
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APPENDIX  A 

COORDINATE  TRANSFORMATIONS 

The  calculations  and  algorithms  presented  in  this  paper  necessitate  describing 
vectors  in  several  different  coordinate  frames.  These  frames  are  all  orthogonal. 
Transformation  of  a  vector  from  one  reference  frame's  coordinates  to  another's  is 
accomplished  using  direction  cosine  matrices  (DCMs). 

Reference  Frames 

The  reference  frames  of  interest  are  as  follows  (see  fig.  A.l ): 

Geocentric  inertial  (IJK)  -  The  I  vector  points  inertially  in  the  vernal  equinox 
direction,  with  I  and  J  lying  in  the  equatorial  plane,  and  K 
piercing  the  north  pole  (coincident  with  earth's  spin  axis). 

Orbit-fixed  (PQW)  -  The  orbital  plane  is  the  fundamental  plane,  with  P 
pointing  to  perigee  (point  of  closest  approach),  Q  is  the  semi- 
latus  rectum  direction,  and  W  is  the  orbit-normal  (also  orbital 
angular  momentum  direction). 

Radar  plane  (XYZ)  -  The  radar  plane  is  the  fundamental  plane,  with  X 
pointing  to  the  maximum  north  latitude  point,  V  lies  in  the 
equatorial  plane,  and  Z  is  the  radar  plane  normal. 

Site-based  (HEN)  -  H  is  local  vertical,  E  is  east,  and  N  is  north.  The  origin 
coincides  with  the  radar  site  location  of  interest. 

Direction  Cosine  Matrices 

Three  DCMs  are  developed  here,  performing  one  orthogonal  rotation  at  a 
time,  then  combining  the  rotations. 
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UK  (Geocentric  inertial) 
K 


PQW  (Orbit-fixed) 


XYZ  (Radar  plane) 


HEN  (Site-based) 


Figure  A.l.  Reference  frames. 

Orbit-fixed  to  Inertial 

Satellite  position  is  most  simply  expressed  in  orbit-fixed  coordinates, 
but  for  comparison  with  earth-fixed  quantities,  it  is  useful  to  transform  both  to  an 
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intermediate  reference  frame.  The  inertial  coordinate  frame  (in  this  case, 
geocentric  inertial)  is  the  obvious  choice.  Fig.  A.2  shows  the  individual  rotations 
required  to  transform  a  vector  from  orbit-fixed  to  inertial  coordinates. 


rotl  = 


cosQ 

sin  Q. 

0 

-sin  Q 

cosQ 

0 

0 

0 

1 

rot  2  = 


1  0  0 

0  cosi        sin  i 

0        -sin  i      cosi 


rot  2  = 


COS  CO 

sin  co 

0 

-sin  co 

cos  CO 

0 

0 

0 

1 

Figure  A.2.  Orbit-fixed  to  inertial  rotations. 
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Combining  the  rotations  yields 

C^  =  [rot3][rot2][rotl] 

ccoc£>  -  scosQci  ccosQ  +  scoc£>ci  scosi 

-scocQ-ccosQci  -scosQ  +  ccocQci  ccosi 

sQsi                       -sicQ  ci 


Because  the  transformation  is  orthogonal, 

c^=[cp^]T  = 


ccocQ  -  scosQci 

ccosQ  +  scocQci 

scosi 


-scocQ-ccosQci      sQsi 
-scosQ  +  ccocQci    -sicQ 
ccosi  ci 


P,     P2     P3~ 
P<     P5     P, 
P,     P8     P9 


(A.l) 


(A.2) 


Inertial  to  Site-based 

Fig.  A. 3  shows  the  individual  rotations  required  to  transform  a 
vector's  coordinates  from  the  inertial  to  a  site-based  reference  frame.  Rotation  #1 
keeps  track  of  earth  rotation.  Rotation  #3  points  out  the  fact  that  the  E  and  N 
vectors  are  not  truly  east-  and  north-pointing.  E  is  in  fact  tangent  to  the  great 
circle  circumscribed  by  the  radar  plane,  and  N  is  orthogonal  to  E  in  the  local 
horizontal  plane. 
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rotl  = 


cos  <J>      sin  <j>     0 

-sin  <J>    cos0    0 

0  0        1 


^l 


where   0=  0)©(t-to)+ lonnle 

(%  is  earth's  rotation  rate 

t0    is  time  of  last  Greenwich  crossing  of  I  axis 

Ion  gllc   is  longitude  of  radar  site 


rot  2  = 


where  J  is  latitude  of  radar  site 


cos  I      0     sin  l 

0         1        0 
-sin  I     0     cos  I 


10  0 

0      cos(az)      sin(az) 
0    -sin(az)    cos(az) 


Figure  A.3.  Inertial  to  site-based  rotations. 
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Combining  the  rotations  yields 
C%  =  [rot3][rot2][rotl] 


clcty  c/s<{)  s/ 

-s0c(az)-s^c(|)s(az)      c<j)c(az)-s(az)s^s<J)      s(az)c/ 
s(|)s(az)-s/c(})c(az)      -c<}>s(az)-s£s<J>c(az)    c(az)c/ 


(A.3) 


Inertial  to  Radar  Plane 

The  transformation  from  inertial  to  radar  plane  coordinates  is  very 
useful  to  the  determination  of  satellite  crossings  of  the  radar  plane.  Rotation  #1 
accounts  for  earth  rotation.  Fig.  A. 4  shows  the  individual  rotations  necessary  to 
perform  this  transformation. 

Combininc  these  rotations  yields 


Ctf  = 


C0     C0 

c0    s0 

-s0 

c0 

0 

-s0   c0    -s0   s0    c0„ 


(A.4) 
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i 

Y 

e 

■^X             rotl  = 

cos  9      sin  6     0 
-sin  6     cos  8     0 
0           0       1_ 

^    1 

where  6  =  co©(t  - 10 )+  lonx 

coe  is  earth's  rotation  rate 

1 0    is  time  of  last  Greenwich  crossing  of  I  axis 

lonx  is  West  longitude  of  maximum  North 

latitude  of  radar  plane 

K 

Z 

) 

* 

% 

^X            rot2  = 

cos6„,      0     sinG_ 

rp                                      rp 

0          1         0 
-sinG^     0    cosG^ 

rp                                    rp 

where  Grp  is  inclination  of  radar  plane  to  equatorial  plane 

Figure  A.4.  Inertial  to  radar  plane  rotations. 
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APPENDIX  B 

LINEARIZED  FILTER  MATRICES  <t>  AND  H 

The  extended  Kalman  filter  can  be  viewed  as  a  predictor-corrector  type 
estimator.  States  are  predicted  using  non-linear  plant  dynamics.  Then  a  correction 
is  computed  by  applying  the  Kalman  gain  to  the  filter  residuals  (observations 
minus  estimates).  Calculation  of  the  Kalman  gain  requires  some  type  of 
linearization  of  the  plant  and  measurement  dynamics  (see  eq.  3.18).  Normally  a 
Taylor  series  expansion  is  used,  keeping  only  first  order  terms  which  are  evaluated 
at  the  best  current  state  estimate.  The  general  form  of  the  linearization  of  a 
nonlinear  discrete  function 


xk+i  =f(xkfTk)  = 


"fi(*k.Tk)' 

f2(xk,Tk) 

fra(xk,Tk) 


(B.l) 


is 


9f,        9f, 


0  = 


3xj       3x 
9f2*        . 

ax~" 


2* 


it 

3x. 


(B.2) 


Bx, 


ax„ 
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Development  of  Linear  Plant  Matrix  O 

As  will  be  seen,  the  state  matrix  O  can  be  as  complex  as  one  lets  it 

become.  Following  is  a  presentation  of  three  progressively  more  complex  models; 

1)  Keplerian  two-body  motion, 

2)  inclusion  of  dominant  first  order  earth  oblateness  gravity  effects,  and 

3)  dominant  atmospheric  drag  effects. 

Two-Body  Problem 

State  dynamics  for  a  Keplerian  2-body  orbit  are  given  as 


Linearizing  gives 


ak+l 

e,.l 

ik+I 
Q  +] 

1  +i 

= 

Wk+I 

Mk+1. 

ek 

_Mk+a;%Tk_ 


=  f(xk,Tk) 


(B.3) 


o  = 


0    0    0    0    0 


0  10    0    0    0 


0  0     10    0    0 


0  0    0     10    0 


0  0    0    0    10 

-a^  Tk     0    0    0    0     1 
2    X-i   k 


(B.4) 
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Secular  Variations  Due  to  Earth 's  Oblateness 

To  first  order,  the  last  three  orbital  elements  are  affected  by  secular 
variations,  so  that 


k+l 


3J 


Qk ^-(cosik)Av 


2p^ 


co,  + 


3J, 


2  —  sin'  i, 


Av 


J 


Mk  +  a' 

k 


l  +  i4^(l-3sin:6k) 

^k 


T, 


=  f(xk,Tk) 


(B.5) 


where 


r,  = 


Pk 


l  +  ek  cos  vk 


sin5k  =  sin  ik  sin(vk  +  cok) 


(B.6) 
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Linearizing  gives 


o  = 


0 
0 


0 


1 


0 

acx 


3ek 

3ek 
dMk+] 


0 


0 


3ik 

3M 


k  +  1 


3e  k  di k 


0 
0 
0 
1 
0 
0 


0 
0 
0 
0 


3M 


k*l 


3co, 


0 


0 


0 

aco^. 
aivi 

9Mk 

"aiv? 


(B.7) 


where  the  terms  in  the  fourth  row  are  found  to  be 


30k+,      3J„  dp 

-f  cosikAv- 


3a, 


Pk 


3a 


aak+;    3J2  _      _aPi 


cosuAv 


aek         pk  3ek 


ank+1    3j2  .  . 

— —  =  — r  sin  i  k  Av 
aik       2p; 


(B.8) 


aok+]    3j2     .   avk 

cost, 


aMk     2Pk2      k  aMk 
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where 


3a  k 

=  l-e! 

9pk 
3ek 

=  -2ae 

3v         l(l +  e2)sinE  +  2esinEcosE 
dMk  -sin  vk(l  -ecosE)J 


The  fifth  row  terms  are  found  to  be 


3C°-      ~3J'A,  Av3Pk 


3ak  pk  3ak 


3W-       "^A.Av^ 


3ek  p\  3ek 

3cok+!      -15J 


2  sin  ik  cosikAv 


3Mk        2p2      J*  3Mk 


where 
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(B.9) 


9ik  2p;  (B.10) 


dcok+1  .    -3J2  A      dvk 


A,    =  2--sin2i  (B.ll) 


Sixth  row  terms  can  be  evaluated  by 

3M_      -9J 


k  +  l 


dL 


3a 


LA'-At-Tkaak 


4  4»        6„       k 


3ek         2rk 
9Mk+1      -9  J 


3ik 


A.    sin6kTk 


3ek 

3(sin5k) 


8i 


dMk+]      -9J2  8(sin5k) 

-  =  — —A6#  sin8kTk- 


8co, 


3co 


3Mk+1  _n  3A    T   3At>    avk 
9Mk     "       2     6°    k   dvk    3Mk 


(B.12) 


where 


3rk 

3a  k 

As 

A7 

3rk 

-2akekA7>  -pk  cosv„ 

3ek 

A27 

8Ag 

-3J2 

4 

3rk      Q   .   2  c     drk 
— — +  3sin    5k  — — 
d\\                      d\\ 

3vk 

3rk 

ak(l-ek  )ek  sin  vk 

3(sin5k) 


dvk  (l  +  ecosvk)' 


d(sin5k) 
3vk 


=  sin  ik  cos(vk  +cok) 


(B.13) 
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and 


A.    =  l-3sin26k 


A,    = 


Y2 
) 


A7    =  l  +  ek  cos \\ 
A*«  =-r(l-3sin28t) 


(B.14) 


Development  of  z  and  H 

The  relationship  between  the  measurements  (pseudo-observables)  and 
orbital  parameters  (filter  states)  is  contained  in  the  nonlinear  function 


zk  =g(x,>TJ 


(B.15) 


and  is  the  same  regardless  of  the  state  dynamics  model  being  considered.    The 
observables  of  interest  here  are  the  north-south  and  east-west  direction  cosines,  or 


z  = 


cos  a 
cosp 


(B.16) 


and  are  related  to  the  range  vector  p  by 


cos  a 
cosp 


Pn. 


(B.17) 
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Note  (from  fig.  B.l)  that  the  HEN  coordinate  system  in  which  p  is  received  is  not 
truly  up,  east  and  north.  H  is  local  vertical,  but  E  is  tangent  to  the  radar  plane 
great  circle,  and  N  is  orthogonal  to  E  in  the  local  horizontal  plane.  This  site- 
dependent  information  is  included  in  the  transformation  matrix  C^1  (see  App.  A) 
as  the  angle  az  between  true  east  and  E. 


Figure  B.l.  Site-based  coordinate  system. 


So  pt,  pN,  and  range  magnitude  p  need  to  be  expressed  in  terms  of  the 
orbital  parameters  in  order  to  form  z  =  g(x,  T).  It  can  be  stated  that 


P  = 


~Ph~ 

pcosy 

Pe 

= 

pcosot 

.Pn. 

pcosP 

(B.18) 
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In  order  to  describe  the  direction  cosines  in  term  of  orbital  parameters,  first 
consider  fig.  B.2.  Earth's  flattening  has  been  exaggerated  in  the  figure  to  clearly 
show  the  effect  on  radar  plane  orientation.  From  the  figure,  it  is  clear  that 


r  =  r 


afftei 


+  R  ,  +p 

Etc         r 


(B.19) 


where 


RsUf  =R(lat)  +  alt, 


R(lat)=    1- 


i-X 


298.26 


298.26 


sin2  (lat) 


(B.20) 


r0ffsC1  -  20km 


Figure  B.2.  Relationship  between  satellite  position  and  range. 
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R(lat)  is  the  latitude  dependent  radius  of  a  point  on  the  reference  geoid  (due  to 
flattening)  [7]  and  is  assumed  to  be  in  the  H  direction.  r0ffseIis  assumed  to  be 
totally  in  the  negative  N  direction,  and  alt^e  is  the  receiver's  altitude.  Combining 
the  elements  of  eqs.  B .  1 8  and  B .  1 9  yields 


rH  = 


R^+pcosy 


pcoscc 


-r 


offset 


+  pcosP 


(B.21) 


rH  is  the  transformation  of  rp  from  the  orbit-fixed  coordinate  system  to  site-based 
coordinates,  or 


H, 

H2 

H3 

P, 

P 

H. 

H5 

H6 

P. 

P 

H7 

Hs 

H9_ 

.P7 

P 

r(R7  cos  v  +  Rg  sin  v) 
r(Rj  cos  v  +  R2  sin  v) 
r(R4  cos  v  +  R5  sin  v) 


3 

|~r  cos  v 

6 

rsin  v 

9  _ 

L  0  . 

(B.22) 


where 


R1=H4P1+H5P4+H6P7 


R2=H4P2+H5P5+H6Pg 


R4  =H7P1+HgP4+H9P7 
R5=H7P2+HgP5+H9Pg 
R7=H1P,+H2P4+H3P7 
Rg=H1P2+H2P5+H3Pg 


(B.23) 
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Combining  eqs.  B.18,B.21,  and  B.22  yields 


r-               -i 

Ph 

Pe 

= 

.Pn. 

r(R7  cos  v  +  Rg  sin  v)-Rale 

r(Rj  cos  v  +  R:  sin  v) 
r(R4  cosv  +  R5  sin  v)  +  r 


oSse; 


(B.24) 


From  eq.  B.24,  p  can  be  calculated  as 

P  =  VPh+Pe+Pn 


(B.25) 


Now,  eq.  B.15  can  be  used  as  the  estimate  of  the  observation  to  compute  filter 

residuals.    But  the  Kalman  gain  depends  on  the  matrix  H,  which  is  derived  by 
linearizing  g(x  k ,  Tk )  with  respect  to  xk. 

The  form  of  H  is  as  follows 


H  = 


afpi 


3a 

a(p-v 


a(P> 


3M 

af,fV 


a.M 


(B.26) 
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where 


3a  p   3a         e  p2  3a 


afp> 


1  3p 


E  1       * 


3M         p  3M 

He  p2  3M 

<VP)    lap. 

3a          p   3a 

1   3p 
p    3a 

<P*0    i*. 

1    3p 

-  n — 

3M         p  3M       '   p2  3M 


Following  are  the  details  of  this  development. 

^Pf      /r.  \  3r 

— -  =  (R ,  cos  v  +  R  2  sin  v)  — 

3a  3a 


dpH      /r.  \  dr 
=  (R7  cos  v  +  Rg  sin  v)  — 

3a  3a 


(B.27) 


3pN ■      ,  x 3r 

— -  =  (R4  cosv  +  R5  sinvj  — 

da  da 


3a     pV       3a  3a  3a 


(B.28) 
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where 


dr         1-e2 
3a      1  +  ecosv 


(B.29) 


The  expressions  for  the  change  in  p  with  respect  to  eccentricity  e  are  identical  to 
eq.  B.28  except  for  the  replacement  of  °*A    with  C^A  ,  which  is 


3r      -a(2e  +  cos  v  +  e 2  cos  v) 
3e  (1  +  ecosv)2 


(B.30) 


Since  the  Rn  terms  in  eq  B.24  are  functions  of  the  orbital  elements  i,  Q,  and  co, 
the  expressions  for  the  change  in  p  with  respect  to  these  three  elements  are  as 
follows,  with  the  replacement  of  3C  by  the  appropriate  element: 


3pE   _fdR1 


cos  V  + 


br, 

3; 


sin  v 


3pF      f3R.  3R„ 

=   cosv  +  — -sinv 

3c      y  3c  3; 


3pN      f3R4  3R5 

cos  v  +  — -sin  v 


3c,      y  3; 


^; 


3a      p^        dq  dq  dq 


(B.31) 
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The  expressions  for  the  change  in  p  with  respect  to  a  change  in  mean  anomaly 

M  is  more  involved  because  of  the  relationship  between  mean  and  true  anomalies. 
This  relationship  is  repeated  here  for  convenience: 


M  =  E  -  e  sin  E 


(B.32) 


COS  V  = 


cosE  -e 
1-ecosE 


(B.33) 


sin  v  = 


inEVl  - 


sin 


1  -ecosE 


(B.34) 


Taking  partial  derivatives  as  appropriate  yields 


dM 


9(cosv)  3(sinv) 

R  7  — z \~  R 


dM 


9M 


r  +  [R7  cos  v  +  R8  sin  v] 


_3r_ 
9M 


dpE 

3M 


'      3(cosv)  |        a(sinv) 


3M 


3M 


r  +  [R,  cos  v  +  R2  sin  v] 


_3r_ 
8M 


3(cosv)  a(sinv) 

K  . r  K  ,  

4      3M  5     3M 


r  +  [R4  cos  v  +  R5  sin  v] 


3M      (B.35) 


ap    i 


5Ph  ,  „    3pE  _    9pN 


aM     pLPh  3M  +Pe  aM+Ps'  3M 
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where 


3(cosv)  .       d\'  dE 

=  -  sin  v 

dM  3E  3M 


3(sin  v)  dv  3E 

cos  v 


dM  3E  9M 

3r  ae(l-e2)sinv  av  3E 

3M  (l  +  ecosv):    3E  8M 

9v  -(l +  e2)sin  E  +  2esinEcosE 
3E  -sin  v(l  -ecosE)2 

3E  _         1 

3M  1-ecosE 


(B.36) 


All  the  expressions  for  the  partial  derivatives  of  the  Rn  terms  can  be  found  in 
Table  B.l. 
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Table  B.l.  DEVELOPMENT  OF  TERMS  FOR  H 


di 

=  scl)s^slH4  -  scoc^siHj 

+  soociH6 

3R2 

ai 

=  ccosQsiH.  -ccocQsiH, 

4                                         5 

+  ccociH 

D 

dR, 

an 

=  -P4H4+P,Hs 

aR: 
an 

=  -P5H4+P2H5 

aR] 

9co 

=  P2H4+P5Hs+PgH6 

aR: 
aw 

=  -P,H4-P4H5-P7H6 

aR4 

9i 

=  scosQsiH-  -scocQsiHs 

+  scociH9 

dR5 

ai 

=  ccosQsiH ,  -  ccocQsiH  B 

+  ccociH, 

dR, 

dQ 

=  -P4H7+P1Hg 

aRi 
an 

=  -P5H7+P2Hs 

3R4 

dco 

=  P2H7+PsHg+P8H9 

3RS 

aw 

=  -P1H7-P4Hg-P7H9 

aR. 

ai 

=  scosQsiH,  -scocQsiH, 

+  scociH3 

ai 

=  ccosQsiHj  -ccocQsiH2 

+  ccociH, 

aR? 
an 

=  -P4H1+P,H2 

dR, 

an 

=  -P5H1+P2H2 

aR7 
aw 

=  P2H,+P5H2+PgH3 

3R8 
aco 

=  -P!H1-P4H2P7H3 
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APPENDIX  C 


FILTER  CODE  (MATLAB) 

Following  is  the  Matlab  code  for  the  extended  Kalman  filter  described  in  this 
paper.  This  particular  filter  models  the  2-body  orbital  dynamics  problem, 
including  secular  and  periodic  perturbations.  It  requires  2  data  files  for  input  (see 

header  of  filter  code): 

datain  -  contains  a  time-tagged  initial  orbital  element  set  and  an  unspecified 
number  of  time-tagged  pairs  of  direction  cosines  and 
corresponding  receiving  sites 

siteinfo  -  contains  pertinent  information  about  the  set  of  receiver  sites 
comprising  the  radar  plane  detection  system. 

The  filter  calls  3  functions: 

enewton  -  performs  a  Newton-Rapson  iteration  to  obtain  eccentric  anomaly, 
given  a  mean  anomaly  and  eccentricity 

intper  -  performs  a  Newton-Rapson  iteration  to  propagate  satellite  position  to 
an  in-radar-plane  condition. 

tH'opi  -  forces  angular  measurements  to  a  number  between  0  and  2n. 

Extended  Kalman  Estimator 

The  code  for  these  functions  follows  the  filter  code. 

This  estimator  requires  2  input   files  to  run 
DATAIN  &  SITEINFO. 

datain=[250   116  82   191626.933  0   0 

103.591322   .00756997   66.81264  189.71773    26.77689 
116.53101 

83    65024.876   6  14.839   .092  0 

83    83436.223   5  -76.912   .012  0 
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83  83436.259  4  -63.554  .009  0 

83  83436.397  1  58.982  .062  0 

83  83436.408  3  -57.058  .043  0 

83  83436.425  6  -78.987  .034  0 

83  160057.722  6  64.179  .030  0 

83  174304.847  3  -16.315  .020  0 

83  174304.904  5  -58.495  .007  0 

83  174304.906  6  -62.516  -.001  0 

83  174304.908  1  74.290  .012   0]/ 

siteinfo=[l  32.57840   -116.97016   1.87473e-5   8.58120 

2  33.44600  -106.99824   2.212449e-4   3.13900 

3  33.33239  -93.55039   8.143107e-6   -4.27563 

4  33.14672  -91.02096   6.738134e-7   -5.66895 

5  32.28768  -83.53628   1.118243e-5   -9.71924 

6  32.04323  -81.92609   3.866064e-6   -10.5800] 

SL*************************************************************** 

function [xkeep, Tkeep, Tint , del Tint , res id, xint , Pkeep] = . . . 

orbper (datain, siteinfo) ; 
inpels= [datain (1,  : )  ' ; datain (2,  : )  ' ]  ; 
observ=datain (3 : length (datain ( : ,  1) )  ,  1 : 5)  ; 
K=inpels (1)  ; 
for  i=l:7 
lat  (i)=siteinfo (i,2) *pi/180; 
Ion (i)=siteinfo (i,  3) *pi/180; 
alt (i)=siteinfo(i,4) ; 
az (i) =siteinfo (i, 5) *pi/180; 
end 

c 

THIS  PROGRAM  PROPAGATES  A  SATELLITE  ORBIT  AS  A  2 -BODY  PROBLEM 
WITH  SECULAR  VARIATIONS  IN  OMEGA,  omega,  and  MEAN  ANOMALY.   IT 
CALLS  THE  FOLLOWING  FUNCTIONS: 

TWOPI  -  places  angle  between  0  &  2*pi. 

INTSEC  -  iterates  to  an  in-fence-plane  condition  when  a  fence 

plane  crossing  is  detected. 
ENEWTON  -  solves  for  eccentric  anomaly,  given  values  for  mean 
anomaly  and  eccentricity. 

*  DEFINE  EARTH  ROTATION  RATE   * 

earthrot=. 05883378171654;   %  Define  earth  rotation  rate  (rad/TU) . 

J2=1.082645e-3; 

a**************************************************************** 

*  DEFINE  VARIOUS  TERMS   * 

£****************************************** 
O 

R=[4e-8  0;0  4e-8];   %  Measurement  noise  matrix. 
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Pkkml=eye (6) *le-8;   %  Initialize  error  covariance  matrix. 

Pkeep ( : , 1 : 6) =Pkkml; 

resid(:,l)=[0;0]; 

G=zeros  (6, 2) ;  %  Initialize  Kalman  gain. 

f inc=33 . 58310*pi/180;     %  Define  fence  inclination. 

cosfinc=cos (fine) ; 

sinfinc=sin (fine) ; 

lonx=-101 . 31348*pi/180;  %  Define  location  of  fence  x-axis. 
<t**************************************************************** 

INITIALIZE  STATES  (ORBITAL  ELEMENTS)    * 
a**************************************************************** 

xnow=[  (inpels (7)  / (2 *pi* 13 . 44  68  9317) ) A (2/3)  ; 

inpels (8) /inpels (9:12) . *pi/180] / 
xkeep  ( : , 1 ) =xnow;   %  Save  original  states. 

S-**************************************************************** 
o 

*  CALCULATE  INITIAL  ECCENTRIC  ANOMALY   * 

o 

EA ( 1 ) =Enewton (xnow (6) , xnow ( 6)  , xnow (2)  )  ; 

S-**************************************************************** 

*  FIND  INITIAL  r  and  TA   * 

0 

cosTA= (cos (EA(1) ) -xnow (2) ) / (1 -xnow (2) *cos  (EA ( i )  )  )  ; 
sir.TA=sqrt  (1-xnow  (2)  "2)  *sin  (EA  (1)  )  /  (1-xnow  (2)  *cos  (EA  (1)  )  )  ; 
r  =  xnow (1) *  (1-xnow  (2) "2) /  (1+xnow (2) *cosTA)  ; 
TA (1 ) =twopi (asin (sinTA) ) ; 
if  TA(1)  <pi 
if  cosTA<0 

TA(l)=pi-TA(l)  / 
end 
else 
if  COSTA<0 

TA(1) =3*pi-TA(l)  / 
end 
end 


%  Perform  quadrant  check 
%     for  true  anomaly. 


a***************************************************************** 
*   CALCULATE  TRANSFORMATION  MATRIX  (ORBITAL  TO  INERTIAL)   * 

o 

cosperi  =  cos (xnow (5)  )  ; 
sinperi=sin (xnow (5)  )  ; 
cosascnd=cos (xnow (4 ) )  ; 
sinascnd=sin (xnow (4  )  )  ; 
costilt=cos (xnow (3) )  ; 
sintilt=sin (xnow (3)  )  ; 


%  Calculate  sines  &  cosines  of 
%  angular  orbital  elements. 


Pl=cosperi*cosascnd-sinperi*sinascnd*costilt ; 
P2=-sinperi*cosascnd-cosperi*sinascnd*costilt / 
P4=cosperi*sinascnd+sinperi*cosascnd*costilt ; 
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P5=-sinperi*sinascnd+cosperi*cosascnd*costilt; 
P  7 =s  inper i  *  s  int  i 1 1 ; 
P8=cosperi*sintilt ; 

S,**x****************************************x********x*********x* 
C 

*  INITIALIZE  TIME   * 

O 

oneJan92=98.9202  5  0  931*pi/180; 
TUperday=l 07. 087 9334 0  97483; 
secperTU=80  6. 81359; 
JD=inpels (3)  ; 
time=inpels (4 )  ; 
hrs  =  fix(time*le-4)  ; 
rnin=fix  (rem  (time,  le4)  *le-2)  ; 
sec  =  rem (time,  le2  )  ; 

Tfilter= (hrs* 3600+min* 60+sec) /secperTU; 
Tgnwch=rem (earthrot*TUperday* (Tf ilter/TUperday+ . . . 

(JD-1) ) +oneJan92, 2*pi) /earthrot; 
Tcld=Tfilter; 
Tkeep  (1) =JD+Tfilter*secperTU/86400;      %  Save  original  time 

*  INITIALIZE  TIME  OE  EIRST  OBSERVATION 

time=observ (1 , 2) ; 

hrs=f ix (time*le-4)  ; 

min=fix (rem (time, le4) *le-2 ) ; 

sec=rem (time, le2 ) ; 

timeobs= (hrs* 3600+min* 60+sec) /secperTU; 

c 

*  CALCULATE  OUT-OE-FENCE-PLANE  COORDINATE  OE  SATELLITE  POSITION  * 

o 

theta=twopi (earthrot *Tgnwch+lonx)  ; 

sintheta=sin (theta) ; 

costheta=cos (theta)  ; 

zfence= .0031+ (-sinfinc*costheta*Pl. . . 

-sinfinc*sintheta*P4  . . . 

+cosf inc*P7 ) *r*cosTA. . .   %  Calculate  current  magnitude 

+  (-sinf inc*costheta*P2 . . . %  of  out-of-plane  coordinate. 

-sinf inc*sintheta*P5 . . . 

+cosfinc*P8) *r*sinTA; 


*   PROPAGATE  ORBITAL  ELEMENTS  TO  NEXT  FENCE-PLANE  INTERSECTION   * 

J2=1.082645e-3; 
p=l; 


81 


n=0; 
for  k=l:K 
Tloop=pi/8*xnow (1) "1 . 5; 

*    INCREMENT  MEAN  ANOMALY  &  APPLY  SECULAR  VARIATION 

mnmotion=xnow (1 ) A (-3/2) ; 

xnow (6) =twopi (xnow (6) +mnmotion*Tloop) ; 

EA  (k  +  1 ) =enewton (xnow ( 6) , xnow ( 6) , xnow (2 ) ) ; 

cosTA= (cos (EA(k+l) ) -xnow (2) ) / (I -xnow (2 ) *cos (EA (k+1 ) ) ) ; 

sinTA=sqrt (1 -xnow (2) "2 ) *sin (EA (k+1 ) ) / (1 -xnow (2) *cos (EA (k+1) ) ) ; 

r=xnow  (1) * (1-xnow (2) "2) / (1  +  xnow (2) *cosTA) ; 

TA (k+1 ) =twopi (asin (sinTA) ) ; 
if  TA(k+l)<pi 
if  cosTA<0 

TA(k+l)  =pi-TA(k  +  l )  /        %  Perfcrir.  quadrant  check 
end  %    for  true  anomaly. 

else 
if  cosTA<0 

TA(k  +  l) =3*pi-TA(k+l)  ; 
end 
end 

£-****xxx*****xxxx****xxxxxx******x***x:**x*x**x 
c 

*    APPLY  SECULAR  VARIATIONS  TO  omega  &  OMEGA    * 

&xxxxxxx*xxxxxxxxxxxxxxxxxxxxxxxxxxxxx**xxxx.xxx 

delTAsec=twopi (TA  (k+1) -TA  (k) ) ; 

semilat=xnow (1 ) *  (1-xnow  (2) ~2) / 

Ql=-3*J2/semilat*2; 

SOMEGA=Ql*cos (xnow (3) ) /2; 

Somega=  (-Q1) *  (2 -5/2* sin (xnow (3) ) A2) /2; 

delsOMEGA=SOMEGA*delTAsec; 

delsomega=Somega*delTAsec; 

xnow (4 ) =xnow (4 ) +delsOMEGA; 

xnow  (5) =xnow (5) +delsomega; 

*  CALCULATE  ATMOSPHERIC  DENSITY   * 

%************************************************xx************** 

k2=3e-6; 
k3=21/6378.135; 
atmdens=k2*exp ( (1-r) /k3)  ; 

*  APPLY  DRAG  PERTURBATIONS  TO  a  &  e   * 

kl=-4*6378.135; 

deldrag= [kl*xnow (1) "2* (1+xnow (2) ) *atmdens*delTAsec; 

kl*xnow (1) * (1-xnow (2) A2) *atmdens*delTAsec; 

0  ;  0 ;  0 ;  0  ]  ; 
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xnow=xnow+deldrag; 

a^xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

*   CALCULATE  TRANSFORMATION  MATRIX  (ORBITAL  TO  INERT I AL)    * 
a**************************************************************** 


cosperi=cos (xnow (5) ) ; 
sinperi=sin (xnow (5) ) ; 
cosascnd=cos (xnow(4)  )  ; 
sinascnd=sin (xnow (4)  )  ; 
costilt  =  cos (xnow (3) )  ; 
sintilt=sin (xnow (3)  )  ; 


%  Calculate  sines  &    cosines  of 
%  angular  orbital  elements. 


£xxXXxxXXXXXXXXXXXXXXXXXXXXXxxxxxxxxxxxxxxxxxxxxxXXXXXXXXXXXXXXXX 
o 

*   UPDATE  TIME   * 

^xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx^xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

Tgnwch=Tgnwch+Tloop; 
Tf i It er=Tfi Iter +T loop; 

o 

*  CALCULATE  OUT-OF-FENCE-PLANE  COORDINATE  OF  SATELLITE  POSITION  * 

&xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx^xxxxxxxxxxxxxxxxxxxxxxx*xxxxxxxx 

theta  (k+1) =twopi (earthrot *Tgnwch+lonx) ; 

sintheta=sin (theta (k+1) ) / 

costheta~cos (theta (k+1) ) / 

zfold=zfence; 

zfence= .0031+ (-sinfinc*costheta*Pl. . . 

-sinf inc*sintheta*P4 . . . 

+cosf inc*P7) *r*cosTA. . .    %  Calculate  current  magnitude 

+  (-sinf inc*costheta*P2 . . .  %  of  out-of-plane  coordinate. 

-sinf inc*sintheta*P5 . . . 

+cosfinc*P8) *r*sinTA; 

Sk**************************************************************** 
o 

*  CHECK  FOR  FENCE  PLANE  INTERSECTION   * 

a-********************************* 

0 

*  NEGATIVE  TO  POSITIVE  ??   * 

%xxxxxxxx*xxxxxxxxxx*xxxxxxxxxxxxxxxxxxxxxxxxx 

if  zfence>0 
if  zfold<0 
n=n+l; 
negtopos=n; 
xnow=xnow' ; 

[Tloop, Tf ilter, Tgnwch, flag, xnow, r, cosTA, . . . 
sinTA, TA(k+l) , theta (k+1 ) , EA (k+1 ) , z fence] 
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inttemp  (xnow, r, cosTA,  .  .  . 

sinTA, TA(k+l) ,  theta (k+1) , EA (k+1) , PI , P2, P4,  ... 
P5, P7, P6, Tgnwch, Tfilter, Tloop) / 
xnow=xnow ' ; 
end 
end 

*  POSITIVE  TO  NEGATIVE  ??   * 
t**************************************************************** 

o 

if  zfence<0 
if  zfold>0 
n=n+l ; 
postoneg=n; 
xnow=xnow '  ; 

[Tloop, Tfilter, Tgnwch, flag, xnow, r, cosTA, . . . 
sinTA, TA(k  +  l)  , theta  (k+1 )  , EA (k  +  1 ) , zfence]  =  .  .  . 
irttemp (xnow, r, cosTA, . . . 

sinTA, TA( k+1) , theta (k+1) , EA (k+1) , PI, P2, P4 , . . . 
P5, P7, P8, Tgnwch, Tfilter, Tloop) / 
xnow=xnow ' / 
end 
end 

TAold=TA(k+l) ; 

%************■*■********************************■**■*******•**■*■******* 

*  SATELLITE  IN  NAVSPASUR  WINDOW  ??   * 

&-****************************■**************<■********************* 
c 

if  flag==l 

9.**************************************************************** 
c 

*  CALCULATE  TIME  FROM  LAST  FILTER  UPDATE  (LAST  OBS)    * 

9-************************************>r*************************** 
0 

delobs=timeobs-Tf ilter; 
while  abs (delobs) <1 
if  length (observ (:, 1 )) >=p 
del obs=timeobs-Tf ilter ; 

delT=delobs; 

delTint (p) =delT*secperTU; 

Tint  (p)  =  (timeobs*secperTU) /3600; 

Tupdate=timeobs-Told; 

if  Tupdate<0 

Tupdate=Tupdate+TUperday; 
end 
Told=timeobs; 
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Tfilter=Tfilter+delT; 
Tgnwch=Tgnwch+delT; 

4»********************************************XX***X*X***XX****** 

*  INCORPORATE  PERIODIC  VARIATIONS   * 

S-x******x**xxx****xx*x**x*X*xx*xx*xx*xxx*xxxxxx**xxxxXxxxxXxxxx*x 

xmean=xnow; 

daper=3* J2*xnow (I) / (2*semilatA2* (1-xnow (2) "2) ) * . . . 

sin  (xnow (3))"2*cos(2*  (xnow (5) +TA (k+1 ) ) ) ; 
deper=3*J2/  (2*semilat "2 )  *(  (1-3/2 *sin  (xnow  (3)  )  A2)  *cosTA+.  .  . 
l/4*sin  (xnow (3) ) A2*cos (2*xnow (5) +TA(k+l)  ) +. . . 
7/ 12* sin (xnow (3) ) ~2*cos (2*xnow(5) +3*TA (k+1) ) ) ; 
diper=-3*J2/ (8*semilat "2 ) *sin (2*xnow (3) ) * . . . 

cos  (2*  (xnow (5) +TA(k+l)  )  )  ; 
dOMper=3* J2/ (4* semilat "2 ) *cos (xnow (3)  )  * .  .  . 

sin  (2* (xnow (5) +TA(k+l) )  ) ; 
dnewper=-3*J2/ (4*semilatA2) * (l-5/2*sin (xnow (3) ) A2) . . . 
*sin  (2*  (xnow (5) +TA(k+l)  )  )  ; 
dtemp=3*J2/ (2*semilat A2) * (l/xnow(2)  *  .  .  . 
(l-3/2*sin  (xnow (3) ) "2) *sinTA- . . . 
l/4*sin  (xnow  (3)  ) "2 -sin (2*xnow(5) +TA (k+1) )  +  . . . 
7/I2*sin  (xnow (3)  ) "2* sin (2*xnow(5) +3*TA (k+1 )))+.. . 
1/2* ( (l-3/2*sin(xnow(3) ) A2) *sin (2*TA(k+l) ) -. . . 
(1-5/2* sin (xnow (3)  )  "2) *sin (2* (xnow (5) +TA  (k+1) ) )  + .  .  . 
3/ 8* sin  (xnow (3) ) "2* sin (2*xnow (5) +4 *TA (k+1 )  )  )  )  / 
domper=dtemp; 
dMAper=dnewper-domper; 

varper= [daper ; deper; diper ; dOMper ; domper; dMAper ] / 
xnow=xnow+varper;  xnow ( 6) =twopi (xnow (6) ) / 
EA (k+1 ) =enewton (xnow ( 6) , xnow (6) , xnow (2) ) / 
cosTA= (cos(EA(k+l)) -xnow (2) ) / (1-xnow (2) *cos (EA (k+1 ) ) ) ; 
sinTA=sqrt (1-xnow (2) ~2) *sin (EA (k+1) ) / (1-xnow (2) *cos (EA (k+1) ) ) ; 
r=xnow (1) * (1-xnow (2) "2) / (l+xnow(2) *cosTA) ; 

TA (k+1) =twopi (asin (sinTA) ) ; 
if  TA(k+l)<pi 
if  cosTA<0 

TA (k+1) =pi-TA (k  +  1)  ;        %  Perform  quadrant  check 
end  %    for  true  anomaly. 

else 
if  cosTA<0 

TA(k+l)=3*pi-TA(k+l)  ; 
end 
end 

*  CALCULATE  ANGULAR  DISPLACEMENT  OF  RECEIVER  SITE  FROM  I -AXIS   * 
a**************************************************************** 
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phi=twopi  (earthrot*Tgnwch+lon  (observ  (p,  3)  )  )  ; 
sinphi=sin (phi) ; 


sinlat  =  sin (lat (observ (p,  3) )  ) 
coslat  =  cos (lat (observ (p,  3)  )  ) 
sinlon=sin (Ion (observ (p, 3) ) ) 
coslon=cos (Ion (observ (p,  3)  )  ) 
cosperi  =  cos (xnow (5)  )  / 
sinperi  =  sin (xnow (5)  )  ; 
cosascnd=cos (xnow (4) ) ; 
sinascnd=sin (xnow (A ) ) / 
costilt=cos (xnow (3) ) ; 
sintilt=sin (xnow (3) ) ; 


Pl  =  cosperi*cosascnd-sir.peri*sinas  end*  cos  tilt; 

P2  =  -sinperi*cosascnd-cosperi  *sinascnd*costilt; 

P4=cosperi*sinascnd+sinperi*cosascnd*costilt; 

P5=-sinperi*sinascnd+cosperi*ccsascnd*costilt ; 

P7=sinperi*sintilt ; 

P8=cosperi*sintilt ; 

sinaz=sin (az (observ (p, 3) ) ) ; 
cosaz=cos (az (observ (p, 3) ) ) ; 

Hl=coslat*cosphi; 
H2=coslat*sinphi; 

H3=sinlat ; 

H4=-sinphi*cosaz-sinlat *cosphi*sinaz; 

H5=cosaz*cosphi-sinaz*sinlat*sinphi; 

H6=sinaz*coslat; 

H7=sinaz*sinphi-cosaz*sinlat*cosphi ; 

H8=-sinaz*cosphi-cosaz*sinlat *sinphi; 

H9=cosaz*coslat; 


R1=H4*P1+H5*P4+H6*P7 
R2=H4*P2+H5*P5+H6*P8 
R4=H7*P1+H8*P4+H9*P7 
R5=H7*P2+H8*P5+H9*P8 
R7=H1*P1+H2*P4+H3*P7 
R8=H1*P2+H2*P5+H3*P8 


a**************************************************************** 
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*  ESTIMATE  MEASUREMENT   * 

£*********************•*********************** 
c 

Rsite  =  sqrt  (1-  (2-1/298 .2  6) /298 . 2  6*sinlat~2) +alt (observ (p, 3) ) ; 
roffset=.0031; 

rhoH=r* (R7*cosTA+R8*sinTA) -Rsite; 

rhoE=r* (Rl*cosTA+R2*sinTA) ; 

rhoN=r* (R4*cosTA+R5*sinTA) +roffset; 

rho  =  sqrt  (rhoHA2+rhoEA2+rhoN~2 )  / 

zest ( : ,p) = [rhoE/rho; rhoN/rho] ; 

size=xnow (1 ) ; shape=xnow (2) ; tilt=xnow (3) ; 
ascnd=xnow (4 ) ;peri=xnow (5) ;MA=xnow (6) ; 

2-**************************************************************** 

c 

*  CALCULATE  COVARIANCE  MATRIX  P   * 

J2=1.082645e-3; 
sinEA=sin (EA (k  +  1)  )  ; 
COSEA=cos (EA (k  +  1)  )  ; 

sindecl=sintilt*sin (TA(k  +  l) +peri)  ; 
semilat=size*  (l-shapeA2)  / 

*   CALCULATE  APPROXIMATE  delta (TRUE  ANOMALY)     * 

mnir1otion=xnow  (1)  A  (-3/2)  ; 

delMA=mnmotion*Tupdate; 

delTA=delMA; 

Alphi  =  -  (l+shape"2) *sinEA+2*shape*sinEA*cosEA; 

A2phi= 1- shape *cosEA / 

A3phi=2-5/2*sintiltA2; 

A4phi=l-3*sindecl"2; 

A5phi=l -shape A2; 

A6phi=sqrt ( l/size+ J2/r A3*A4phi ) ; 

A7phi=l+shape*cosTA; 

dOMEGAda=3*J2/ (size*semilat~2) *costilt*delTA; 

dOMEGAde=-6*J2*size*shape/semilatA3*costilt*delTA; 

dOMEGAdi=3*J2/ (2*semilat "2 ) *sintilt*delTA; 

dOMEGAdM=-3*J2/ (2*semilat "2) *costilt*Alphi/ (sinTA*A2phiA3) ; 

domegada=-3*J2/ (size*semilat A2) *A3phi*delTA; 
domegade=6* J2* size* shape /semilat A3*A3phi *delTA; 
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domegadi=-15*J2/ (2*semilat ~2) *sintilt*costilt*delTA; 
domegadM=-3*J2/ (2*semilat"2) *A3phi*Alphi/ (sinTA*A2phiA3) / 

drkdak=A5phi /A7phi ; 
dMAda=3/2*Tupdate*A6phi* (-1 /size"2+ J2*A4phi* (-3/rA4) *drkdak) / 

drkdek= (-2*size*shape* (l-shape*cosTA) -semilat *cosTA) /A7phi~2; 
dMAde=3/2*Tupdate*A6phi*J2*A4phi* (-3/r~4) *drkdek; 

dsindkdi=sin (TA(k+l) +peri) * cost i It; 
dMAdi=3*Tupdate*A6phi* (-3*J2/rA3) *sindecl*dsindkdi; 

dsindkdo=sintilt*cos (TA(k+l) +peri)  ; 
dMAdomega=3*Tupdate*A6phi* (-3*J2/r~3) *sindecl*dsindkdo; 

drdTA= semilat* shape* sinTA/A7phiA2; 

dsinddTA=dsindkdo; 

dmessdTA=-3*J2/r~2*drdTA- (-9* J2/r "4 *sindecl"2*drdTA+ . . . 
6*J2/r"3*sindecl*dsinddTA) ; 

dTAdMA=-Alphi/ (sinTA*A2phi"3) ; 

dMAdM=l+3/2*Tupdate*A6phi*dmessdTA*dTAdMA; 

dada=l+2*kl*size* (1+shape) *atmdens*delTA; 
dade=kl*size~2*atmdens*delTA; 
dadM=-kl*size~2* (1+shape) *atmdens*dTAdMA; 

deda=kl*A5phi*atmdens*delTA; 
dede-1-2 *kl *  size* shape *atmdens* delTA; 
dedM=-kl*size*A5phi*atmdens*dTAdMA; 

phimat=[dada  dade  0  0  0  dadM; 
deda  dede  0  0  0  dedK; 
0  0  10  0  0; 

dOMEGAda  dOMEGAde  dOMEGAdi  1  0  dOMEGAdM; 
domegada  domegade  domegadi  0  1  domegadM; 
dMAda  dMAde  dMAdi  0  dMAdomega  dMAdM] ; 

9-******************************************* 
o 

*   CALCULATE  PLANT  NOISE   * 

J2=1.0826e-3; 

semilat=xnow (1) * (1-xnow (2) A2)  ; 

siga2=  (3*J2*xnow (1) / (2* semilat "2* (1-xnow (2)  "2)  )  *  .  .  . 

sin (xnow(3) ) "2* cos (2* (xnow(5) +TA(k+l) ) ) ) "2; 
sige2= (3*J2/ (2* semilat "2) * ( (l-3/2*sin (xnow(3) ) "2) *cosTA+. . . 
sigi2= (3*J2/ (8*semilat A2) *sin (2*xnow(3) ) *. . . 

cos (2* (xnow(5)+TA(k+l) ) ) ) "2; 
mnmotion=xnow  (1) A  (-3/2) ; 
MApert=mnmotion*3*J2/ (2* semilat ~2) *(l-3/2*sin (xnow(3) ) "2) * . . . 

sqrt (1-xnow (2)  "2)  ; 
sigOMEG2= (3* J2/ (2*semilat ~2) *mnmotion*cos (xnow (3) ) *Tupdate) "2; 
sigomeg2= (3*J2/ (2*semilatA2) *mnmotion* . . . 
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(2-5/2*sin  (xnow (3) ) "2) *Tupdate) "2/ 
sigMA2= (MApert*Tupdate) "2; 
Q=[le-2*siga2  0  0  0  0  0 
0  le-3*sige2  0  0  0  0 
0  0  sigi2  0  0  0 
0  0  0  sigOMEG2  0  0 
0  0  0  0  sigomeg2  0 
0  0  0  0  0  sigMA2] .*le-3; 

o 

*  PROPAGATE  ERROR  COVARIANCE   * 

o 

Pkk=Pkkml; 

Pkkml=phimat*Pkk*phimat ' +Q; 

o 

*   CALL  IN  OBSERVATION  DATA    * 

cosalpha=sin (observ (p, 4) *pi/180) ; 
cosbeta=sin (observ (p,5)*pi/180); 
z ( :  ,  p)  =  [cosalpha; cosbeta] / 

9***************************************************************** 
c 

*  CALCULATE  LINEARIZED  MEASUREMENT  MATRIX  H   * 

a.******************************************** 

al=l-xnow(2) A2; 
a2=l+xnow (2) *cosTA; 
a3=2*xnow(2) +cosTA* (l+xnow(2) A2) ; 
a4=2*xnow (2) *sinEA*cosEA; 

a5=l+xnow (2) A2; 
a6=-a5*sinEA+a4  ; 
a7=l-xnow (2) *cosEA; 

drda=al/a2; 
drde=-xnow (1) *a3/a2A2; 

dcosTAdT=-sinTA; 

dTAdEA= (-a5*sinEA+a4)  /  (-sinTA*a7 "2 )  ; 

dEAdMA=l/a7; 

dcosTAdM=dcosTAdT*dTAdEA*dEAdMA; 

dsinTAdT=cosTA; 
dsinTAdM=dsinTAdT*dTAdEA*dEAdMA; 

drdTA=xnow(l) * xnow (2) *al*sinTA/a2 A2; 
drdMA=drdTA*dTAdEA*dEAdMA; 

dRldi=sinperi*sinascnd*sintilt*H4- . . . 

sinperi*cosascnd*sintilt*H5+ . . . 
sinperi*costilt*H6; 
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dR2di=cosperi*sinascnd* 

cosperi*cosascn 

ccsperi*costilt 

dR4di=sinperi*sinascnd* 

sinperi*cosas 

sinperi*costi 

dR5di=cosperi*sinas 

cosperi*cosas 

cosperi*costi 

dR7di=sinperi*sinas 

sinperi*cosas 

sinperi*costi 

dR8di=cosperi*sinas 

cosperi*cosas 

cosperi*costi 


sinti 

d  *  s  i  n 

*H6; 

sinti 

cnd*s 

lt*H9 

cnd*s 

cnd*s 

lt*H9 

cnd*s 

cnd*s 

lt*H3 

cnd*s 

cnd*s 

lt*H3 


lt*H4- . . . 
tilt*H5+. . 

lt*K7-. . . 
intilt*H6+ 

r 

intilt*H7- 

intilt*H8+ 

intilt*Hl- 

intilt*H2+ 

intilt*Hl- 

intilt*H2+ 


dRldOM=-P4*H4+Pl*H5 
dR2dOM=-P5*H4+P2*H5 
dR4dOM=-P4*H7+Pl*H8 
dR5dOM=-P5*H7+P2*H8 
dR7dOM=-P4*Hl+Pl*K2 
dR8dOM=-P5*Hl+P2*K2 


dRldom=P2*H4+P5*H5+P8*H6; 
dR2dom=-Pl*H4-P4*H5-P7*H6 
dR4dom=P2*H7+P5*H8+P8*H9; 
dR5dom=-Pl*H7-P4*H8-P7*H9 
dR7dom=P2*Hl+P5*H2+P8*H3; 
dR8dom=-Pl*Hl-P4*H2-P7*H3 


dpHda= (R7*cosTA+R8*sinTA) *drda; 
dpEda= (Rl*cosTA+R2*sinTA) *drda; 
dpNda= (R4*cosTA+R5*sinTA) *drda; 
dpda= (rhoH*dpHda+rhoE*dpEda+rhoN*dpNda) /rho; 

dpHde= (R7*cosTA+R8*sinTA) *drde; 
dpEde= (Rl*cosTA+R2*sinTA) *drde; 
dpNde= (R4*cosTA+R5*sinTA) *drde; 
dpde= (rhoH*dpHde+rhoE*dpEde+rhoN*dpNde) /rho; 

dpHdi= (dR7di*cosTA+dR8di*sinTA) *r; 
dpEdi= (dRldi*cosTA+dR2di*sinTA) *r; 
dpNdi= (dR4di*cosTA+dR5di*sinTA) *r; 
dpdi= (rhoH*dpHdi+rhoE*dpEdi+rhoN*dpNdi) /rho; 

dpHdOM= (dR7dOM*cosTA+dR8dOM*sinTA) *r; 
dpEdOM= (dRldOM*cosTA+dR2dOM*sinTA) *r; 
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dpNdOM= (dR4dOM*ccsTA+dR5dOM*sinTA) *r; 

dpdOM= (rhoH*dpHdOM+rhoE*dpEdOM+rhoN*dpNdOM) /rho; 

dpHdom= (dR7dom*cosTA+dR8dom*sinTA) *r; 
dpEdom= (dRldom*cosTA+dR2dom*sinTA) *r; 
dpNdom= (dR4dom*cosTA+dR5dom*sinTA) *r; 
dpdom=  (rhoH*dpHdom+rhoE*dpEdom+rhoN*dpNdom)  /rho; 

dpHdMA= (R7*dcosTAdM+R8*dsinTAdM) *r+ (R7*cosTA+R8* 

sinTA) *drdMA; 

dpEdMA= (Rl*dcosTAdM+R2*dsinTAdM) *r  + (Rl*cosTA+R2* 

sinTA) *drdMA; 

dpNdMA= (R4*dcosTAdM+R5*dsinTAdM) *r+ (R4*ccsTA+R5* 

sinTA) *drdMA; 

dpdMA= (rhoH*dpHdMA+rhoE*dpEdMA+rhoN*dpNdMA) /rho; 

drhoEda=dpEda/rho-rhoE*dpda/rho/v2; 

drhoNda=dpNda/rho-rhoN*dpda/rho/v2; 

drhoEde=dpEde/rho-rhoE*dpde/rhCv2; 
drhoNde=dpNde/rho-rhoN*dpde/rho~2  ; 

drhoEdi=dpEdi/rho-rhoE*dpdi/rho"2; 
drhoNdi=dpNdi/rho-rhoN*dpdi/rhoA2; 

drhoEdOM=dpEdOM/rho-rhoE*dpdOM/rho~2; 
drhoNdOM=dpNdOM/rho-rhoN*dpdOM/rhoA2; 

drhoEdom=dpEdom/rho-rhoE*dpdom/rho"2; 
drhoNdom=dpNdom/rho-rhoN*dpdom/rho'N2; 

drhoEdMA=dpEdMA/rho-rhoE*dpdMA/rho"2; 
drhoNdMA=dpNdMA/rho-rhoN*dpdMA/rho"2; 

H=[drhoEda  drhoEde  drhoEdi  drhoEdOM  drhoEdom  drhoEdMA; 

drhoNda  drhoNde  drhoNdi  drhoNdOM  drhoNdom  drhoNdMA] ; 
^•*************************************************** 

*   CALCULATE  KALMAN  GAIN   * 

G=Pkkml*H' *inv(H*Pkkml*H'+R) ; 
res id ( : ,p)=z (:,p) -zest (:,p)  ; 


*   CALCULATE  IMPROVED  COVARIANCE   * 
Pkk=(eye(6) -G*H) *Pkkml; 
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Pkeep  (: ,p*6+l:6* (p+1)  ) =Pkk; 

C^XXXXXXXXXXXXX-XXXXX^RXXXXXXXXX^XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

c 

CALCULATE  IMPROVED  MEAN  ORBITAL  ELEMENTS 

xnow=xmean; 
G* res id ( : ,  p) ; 

xnow=xnow+G*resid ( : ,p) ; 
xnow (6) =twopi (xnow (6) )  ; 

xint (1 : 2, p) =xnow (1:2) ; 

xint (3:6,p) =xnow(3: 6) *  180 /pi; 

*   CALCULATE  IMPROVED  ECCENTRIC  ANOMALY   * 

S-*xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
o 

EA  (k-t-1 )  =enewton  (xnow  (6)  ,  xnow  ( 6)  ,  xnow  (2  )  )  ; 

G 

SOLVE  EOR  IMPROVED  cos (TA)  ,  sin (TA)  ,  r,  sin(decl),  n  * 

^XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

o 

cosTA= (cos(EA(k+l)) -xnow (2) ) / (1-xnow (2) *cos (EA (k+1) ) ) ; 

sinTA=sqrt (1-xnow (2) A2) *sin (EA (k+1) ) / (1-xnow (2) . . . 

*cos (EA(k+l) ) ) / 
r  =  xnow  (1 )  x  (1-xnow  (2) *2)  /  (1+xnow  (2)  *cosTA)  / 

TA (k+1) =twopi (asin (sinTA) ) / 

if  TA(k+l)<pi 

if  cosTA<0 
TA (k+1) =pi-TA (k+1) ;        %  Perform  quadrant  check 

end  %    for  true  anomaly. 

else 

if  cosTA<0 
TA(k  +  l) =3*pi-TA(k  +  l)  ; 

end 
end 

TAold=TA(k+l) ; 


9-xxxxxxxxxxxxxxxxxxxx*xx**xxxxxxxxxxxxxxxxxxxxxx 
o 

if  length (observ (:,  1 ))  >p 
P=P+1; 

time=observ (p, 2 ) ; 
hrs=fix (time*le-4) ; 
min=f ix (rem (time, le4) *le-2) ; 
sec=rem (time, le2) ; 

timeobs=  (hrs*3600+min*60+sec) /secperTU; 

delobs=timeobs-Tf ilter; 
else 

delobs=2; 
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end 
end 
end 
end 
flag=0; 

a**************************************************************** 
*   RESET  TIMES  FOR  NEW  DAY   * 

if  Tgnwch>106. 7955369153698 

Tgnwch=Tgnwch-10  6.7  9553  69153698; 
end 

if  Tfilter>107. 0679334 
Tfilter=Tfilter-107.0  87  9334; 
JD=JD+1; 
end 
xkeep ( : , k+1) =xnow; 

Tkeep (k+1) = JD+Tf ilter*secperTU/8  64 00/ 
end 


In-plane  Iteration 

function [T, Tfilter, Tgnwch, flag, xnow, r, cosTA, . . . 

sinTA, TA, theta, Enow, z fence] = . . . 

inttemp (xnow, r, cosTA, sinTA, TA, theta, Enow, . . . 

Pl,P2,P4,P5,P7,P8,Tgnwch, Tfilter, T) ; 
xnow=xnow ' ; 
J2=1.082645e-3; 

earthrot=. 05883378171 654; %Define  earth  rotation  rate (rad/TU) . 
lonx=-101 .  3134  8*pi/180; 
finc=33.58310*pi/180; 
sinfinc  =  sin (fine)  / 
cosfinc  =  cos (fine)  ; 
delT=0; 
zfence=l;   %  INITIALIZE  OUT-OF-FENCE-PLANE  VALUE  OF  POSITION 

&*****•**************************************** 
O 

%       *    ITERATE  TO  IN-PLANE  CONDITION   * 

while  abs (zfence) >=le-8 
%         *    ADD  IN  PERIODIC  VARIATIONS    * 

o 

semilat=xnow (1) * (1-xnow (2) ~2)  ; 

daper=3*J2*xnow(l) / (2*semilat"2* (1-xnow (2)  "2)  )  * . . . 
sin (xnow(3) ) ~2*cos (2*  (xnow(5) +  TA) ) ; 
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deper=3*J2/  (2* serr.il at "2 )  *  (  (1-3/2* sin  (xnow(3)  )  "2)  *cosTA+.  .  . 
l/4*sin (xnow (3) ) "2* cos (2*xnow (5) +  TA) +  .  .  . 
7/ 12* sin  (xnow  (3)  ) "2* cos  (2*xnow(5) +3*TA)  )  ; 
diper=-3*J2/ (8*semilat ~2) *sin (2*xnow(3) ) . . . 
*cos (2* (xnow(5) +  TA) ) / 
dOMper=3* J2/ (4 *semilat "2 ) *cos (xnow (3) ) *sin (2* (xnow (5) +  TA) ) ; 
dnewper=-3*J2/ (4  *  semi 1 at "2) *(l-5/2*sin  (xnow  (3) ) "2)  . . . 

*sin  (2*  (xnow (5) +  TA)  )  ; 
dtemp=3*J2/ (2*semilatA2) * (1 /xnow (2) *. . . 
( (1-3/2 *sin (xnow (3) ) "2) *sinTA-. . . 
1/4* sin (xnow (3) ) "2* sin (2*xnow(5) +TA)+. . . 
7/ 12* sin  (xnow  (3)  )  "2* sin  (2*xnow(5)  +3*TA)  )  +  .'.  . 
l/2*((l-3/2*sin (xnow (3) )*2)*sin(2*TA)-... 
(1-5/2* sin (xnow (3) ) "2) *sin (2* (xnow (5)+TA) ) +. . . 
3/ 8* sin (xnow  (3)  ) "2* sin (2*xnow(5) +4*TA)  ) ) ; 
domper=dtemp; 
dMAper=dnewper-domper ; 

varper= [daper ; deper ; diper ; dOMper ; domper ; dMAper ] ; 
xper=xnow+varper; 

EAper=enewton (xper ( 6) , xper ( 6)  , xper (2 )  )  ; 

ccsTAper= (cos (EAper) -xper (2) ) / (1-xper (2) *cos (EAper) ) ; 

sinTAper=sqrt (1-xper (2) "2) *sin (EAper) / (1-xper (2) . . . 

*cos (EAper) ) / 
rper=xper  (1) *  (1-xper (2) "2) / (1  +  xper (2) *cosTAper)  ; 

TAper=twopi (asin (sinTAper) )  ; 
if  TAper<pi 

if  cosTAper<0 

TAper=pi-TAper ;        %  Perform  quadrant  check 
end  %    for  true  anomaly. 

else 

if  cosTAper<0 

TApe  r  =  3 *  p i - TApe  r ; 
end 
end 

9-*************************************************** 
o 

RECALCULATE  PERIODIC  TRANSFORMATION  MATRIX    * 
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P7per  =  sin (xper (5)  ) *sin  (xper (3) )  ; 
F8per=cos (xper (5) ) *sin (xper (3)  )  ; 
stuf f=sinfinc*earthrot*xper (1) "1 . 5*rper* . . . 

cosTAper* (Plper*sin  (theta)  . . . 

-P4per*cos (theta) ) +sinfinc*earthrot *xper  ( 1 ) "1.5*.  .  . 

rper*sinTAper* (P2per*sin (theta) -P5per* . . . 

cos  (theta) ) +xper (1) / (1-xper (2) *cos (EAper) ) * . . . 

(-sin (EAper) * (-sinf inc*Plper*cos (theta) - . . . 

sinfinc*P4per*sin  (theta)  +cosfinc*P7per)+.  .  . 

sqrt (1-xper (2) A2) *cos (EAper) * (-sinfinc*P2per* . . . 
cos  (theta)  -sinf  inc*P5per*sin  (theta)  +  cosf  ir.~xP8per )  )  ; 

zfence= .0031+ (-sinfinc*cos (theta) *Plper-sinfinc* . . . 

sin (theta) *P4per+cosf inc*P7per ) *rper* .  .  . 

cosTAper+ (-sinf inc*cos (theta) *P2per. . . 

-sinf inc*sin (theta) *P5per+cosf inc*P8per ) * . . . 

rper*sinTAper ; 
delMA=-z fence/ stuff; 
mnmotion  =  xnow (1 ) * (-3/2)  / 
del2T=delMA/mnmotion; 

c 

PROPAGATE  MEAN  ORBITAL  ELEMENTS  TO  ITERATIVE  VALUE  OF  TIME   * 

2-********************************* 
c 

xnow  (  6) =twopi (xnow (6) +delMA) ; 

Enow=enewton (xnow ( 6) , xnow ( 6)  ,  xnow (2  )  )  ; 

TAold=TA; 

cosTA= (cos (Enow) -xnow (2)  ) /  (1 -xnow (2)  *cos  (Enow)  )  ; 

sinTA=sqrt (l-xnow(2) "2) *sin (Enow) / (1-xnow (2) . . . 

*cos  (Enow)  ) ; 
r=xnow  (1) *  (1-xnow (2) "2) / (1+xnow (2) *cosTA) / 

TA=twopi (asin (sinTA) ) ; 
if  TA<pi 

if  cosTA<0 

TA=pi-TA;        %  Perform  quadrant  check 
end  %    for  true  anomaly. 

else 

if  cosTA<0 

TA=3*pi-TA; 
end 
end 

%     *    APPLY  SECULAR  VARIATIONS  TO  omega  &  OMEGA    * 
%*************************•********************************** 

delTA=TA-TAold; 
if  delTA>pi/2 

delTA=delTA-2*pi; 
end 
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if  delTA<-pi/2 
delTA=delTA+2*pi; 

end 
semilat=xnow (1) * (1-xnow (2) ~2) ; 
Ql=-3*J2/semilat"2; 


SOMEGA=Ql*cos (xnow (3)  )  /2; 

Somega=-Ql* (2 -5/2* sin (xnow (3) ) "2) /2; 

delsOMEGA=SOMEGA*delTA; 

de Is omega =S omega* del TA; 

xnow (4) =xnow (4) +delsOMEGA; 

xnow (5) =xnow (5) + dels omega; 

*  CALCULATE  ATMOSPHERIC  DENSITY   * 

k2=3e-6; 
k3=21/6378.135; 
atmdens=k2*exp ( (1-r) /k3) ; 

*  APPLY  DRAG  PERTURBATIONS  TO  a  S  e   * 

c 

kl=-4*6378.135; 

deldrag= [kl*xnow  (1) ^2* (1+xnow (2) ) *atmdens*delTAsec; 

kl*xnow (1) *  (1-xnow (2) "2) *atmdens*delTAsec; 

0  ;  0 ;  0  ;  0  ]  ; 
xnow=xnow+deldrag; 

o 

*  RECALCULATE  TRANSFORMATION  MATRIX  (ORBITAL  TO  INERTIA!) 

Pl=cos (xnow (5) ) *cos (xnow (4 ) ) - . . . 

sin (xnow (5) ) *  sin (xnow (4) ) *cos  (xnow  (3)  ) 
P2=-sin (xnow (5) ) *cos (xnow (4 ) ) - . . . 

cos (xnow  (5)  ) *sin  (xnow (4 )  ) *cos  (xnow (3)  ) 
P4=cos  (xnow (5) ) *  sin (xnow (4) ) + . . . 

sin  (xnow  (5) ) *cos (xnow (4 )  ) *cos  (xnow  (3) ) 
P 5  = -sin  (xnow (5) ) *  sin (xnow (4) ) +. .  . 

cos  (xnow  (5) ) *cos (xnow (4 ) ) *cos (xnow (3) ) 
P7=sin (xnow (5) ) *sin (xnow (3) ) ; 
P8  =  cos (xnow (5) ) *  sin (xnow (3) ) ; 

o 

%      *      IMPLEMENT  TIME  CORRECTION     * 

delT=delT+del2T; 

theta=twopi (earthrot* (Tgnwch+delT) +lonx) ; 
end 

O 

%       *      UPDATE  ACTUAL  TIME  OF  INTERSECTION     * 
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Tgnwch=Tgnwch+delT; 
Tfilter=Tfilter+delT; 

T-T+delT; 


S.xx*XxxxxxxX*xxxXxxxxxxxxxxxxXxxxxx*xxx*xxxxxXxxXxxxxxxxxxxxx 
o 

%       *    CALCULATE  FENCE-PLANE  COORDINATES     * 

9-************************************************* 
o 

x=  (cosfinc*cos  (theta) *Pl+cosfinc*sin (theta) *P4+  . . . 

sinfinc*P7) *r*cosTA+ (cosfinc*cos (theta) *F2+. . . 

cosfinc*sin (theta) *P5+sinf inc*P8) *r*sinTA; 
y=  (-sin  (theta) *Pl  +  cos (theta) *P4) *r*cosTA. . . 

+ (-sin (theta) *P2+cos (theta) *P5) *r*sinTA; 

<ixxxxxx*xxxxxxxXxxx*xxxx**^*******x^*xxxxxxxxxxxxxxxxxxxxxxxx 

DETERMINE  WHETHER  ANY  RCVR  MAY  BE  ABLE  TO  OBSERVE  SATELLITE  * 

O 

if  x/sqrt  (xA2+yA2) >=cos (40*pi/180) 

flag=l  ; 

end 
zfence=0; 


Eccentric  Anomaly  Iteration 

function  [EA] =Enewton (EA, MA, e) 

if  MA<EA 

if  EA>1 .5*pi 

MA=MA+2*pi; 
end 
end 
errE=l; 

while  abs  (errE) >=le-10 
errE=EA-e*sin (EA) -MA; 
delEA= (-1/ (l-e*cos (EA) ) ) *errE; 
EA=twopi (EA+delEA) ; 
MA=twopi (MA) ; 
end 
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Angle  Between  0  &  2n 


%  This  function  will  take  any  angle  (in  radians)  and 

calculate 

%  its  equivalent  between  zero  and  2*pi 

function  phi=twopi (phi) 

if  phi>=2*pi 

while  phi>=2*pi 
phi=phi-2*pi; 
end 
else 

while  phi<0 

phi=phi+2*pi ; 
end 
end 
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